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Abstract 

Surface  effects  limit  the  performance  of  hypersonic  vehicles,  micro-electro- mechan¬ 
ical  devices,  and  directed  energy  systems.  This  research  develops  methods  to  predict 
adsorption,  scattering,  and  thermal  desorption  of  molecules  on  a  surface.  These  meth¬ 
ods  apply  to  physisorptive  (adsorption  and  scattering)  and  chemisorptive  (thermal 
desorption)  gas-surface  systems.  Engineering  and  design  applications  will  benefit 
from  these  methods,  hence  they  are  developed  under  the  Direct  Simulation  Monte 
Carlo  construct. 

The  novel  adsorption  and  scattering  contribution,  the  Modified  Kisliuk  with  Scat¬ 
tering  method,  predicts  angular  and  energy  distributions,  and  adsorption  probabil¬ 
ities.  These  results  agree  more  closely  with  experiment  than  the  state-of-the-art 
Cercignani-Lampis-Lord  scattering  kernel.  Super-elastic  scattering  is  predicted.  Gas- 
adlayer  interactions  are  included  for  the  first  time.  Accommodation  coefficents  can 
be  determined  by  fitting  simulations  to  experimental  data. 

The  new  thermal  desorption  model  accurately  calculates  angular,  translational, 
rotational,  and  vibrational  distributions,  and  the  rotational  alignment  parameter. 
The  model  is  validated  by  comparing  with  experiments.  Multiple  transition  states 
are  considered  in  a  set  of  non-dimensionalized  equations  of  motion,  linked  with 
temporally-accurate  event  timing.  Initial  conditions  are  chosen  from  a  new  trun¬ 
cated  Maxwell- Boltzmann  distribution.  Run  times  are  improved  by  eliminating  the 
Gaussian  Weighting  of  desorbing  products.  The  absorption  energy  barrier  is  shown 
to  significantly  contribute  only  to  the  translational  energy  of  desorbing  molecules  by 
contributing  energy  to  each  adatom  in  a  similar  manner. 
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Scattering,  Adsorption,  and 
Langmuir-Hinshelwood  Desorption  Models 

FOR  PHYSISORPTIVE  AND  CHEMISORPTIVE 

Gas-Surface  Systems 

I.  Introduction 

Pulse  shortening,  a  negative  effect  which  limits  the  energy  output  of  High-Power 
Microwave  (HPM)1  devices,  has  been  identified  as  the  most  important  problem 
in  the  HPM  held  [16].  Pulse  shortening  occurs  when  the  pulse  generated  by  the  HPM 
device  is  shorter  than  the  design  allows,  effectively  limiting  the  amount  of  energy  that 
can  be  produced  at  the  same  power  level.  The  main  contributors  to  pulse  shortening 
are  plasma  effects,  one  of  which  is  gap  closure.  Gap  closure  is  the  process  by  which 
the  surface  plasma  migrates  into  the  cavity,  thereby  changing  the  diode  impedance. 
This  change  in  turn  reduces  wave  coupling  in  the  Slow  Wave  Structure  (SWS),  a 
critical  step  in  creating  microwaves.  The  plasma  created  at  the  anode  surface  is  a 
result  of  gas  desorption  [both  thermal  desorption  and  Electron-Stimulated  Desorption 
(ESD)],  held  emission,  and  Secondary  Electron  Emission  (SEE).  Being  able  to  predict, 
minimize,  or  prevent  gas  desorption  altogether  thus  becomes  an  important  aspect  in 
advancing  the  state-of-the-art  in  HPM  technology. 

Gas  desorption  is  also  an  important  process  in  other  fields  where  surface  effects 
dominate,  such  as  catalysis,  thin-Elm  processing,  space  environments,  and  hyper- 
sonics.  Unfortunately,  desorption  is  complex  and  difficult  to  study  experimentally, 
theoretically,  and  computationally.  The  many-body  problem  of  desorption  has  no 
^or  a  full  listing  of  all  the  acronyms  used  in  this  document,  please  refer  to  Appendix  B. 
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straightforward  analytical  solution.  Experimentalists  continue  to  seek  methods  of 
measuring  surface  kinetics  (e.g.  orientations,  internal  and  kinetic  energies,  and  elec¬ 
tronic  states).  Even  quantum-based  calculations,  such  as  Density  Functional  Theory 
(DFT),  must  truncate  the  number  of  molecular  interactions  to  achieve  reasonable 
results.2 

Researchers  at  the  Air  Force  Research  Laboratory  (AFRL)  develop  state-of-the-art 
ffPM  systems  and  technologies  through  experimentation  and  simulation.  A  success¬ 
ful  tool,  the  Improved  Concurrent  Electromagnetic  Particle-In-Cell  (ICEPIC)  code, 
allows  researchers  to  both  understand  more  explicitly  what  is  occurring  within  the 
entire  system,  and  quickly  assess  new  ffPM  system-level  designs,  saving  valuable 
resources. 

In  order  to  properly  simulate  the  performance  of  an  HPM  system,  it  is  important 
to  model  the  phenomena  occuring  within  the  magnetron.  For  HPM  systems,  EM 
waves  are  generated  by  a  source,  guided  by  a  waveguide,  and  then  broadcast  (or 
directed)  with  an  antenna.  A  magnetron  is  a  typical  microwave  source.  Microwaves 
are  generated  in  a  magnetron  from  the  interaction  of  moving  electrons  with  both  an 
applied  magnetic  field  and  side  cavities. 

AFRL  is  seeking  to  employ  Direct  Simulation  Monte  Carlo  (DSMC)  techniques 
in  ICEPIC  to  model  both  gas-surface  effects  and  gas-gas  collisions,  including  pulse 
shortening  due  to  surface  contaminants.  DSMC  is  a  popular  engineering  method 
that  approximates  the  solution  to  the  Boltzmann  equation.  The  Boltzmann  equa¬ 
tion  describes  the  position  and  velocity  probability  distributions  of  gas  molecules  in 
both  space  and  time.  Typically,  the  Boltzmann  equation  is  applied  to  flows  with  low 
densities  to  keep  computational  costs  low.  Thus,  DSMC  is  useful  in  the  fields  of  hyper- 
sonics,  Micro- Electro-Mechanical  Systems  (MEMS),  and  vacuum  devices  (e.g.  mag- 

2  Please  note  that  here  and  in  the  remainder  of  this  paper  the  term  molecule  refers  to  an  atom  or 
to  a  molecule. 
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netrons),  to  name  a  few.  DSMC  approximates  the  solution  to  the  Boltzmann  equation 
by  simulating  molecule-molecule  interactions  on  a  coarse  level.  A  number  of  simu¬ 
lated  molecules  are  followed  as  they  interact,  and  through  well-established  DSMC 
techniques,  excellent  solutions  can  be  obtained.  These  collisions  are  between  sim¬ 
ulated  molecules,  each  of  which  represent  many  real  molecules.  Since  system-level 
simulations  are  performed  by  ICEPIC,  an  engineering  method  such  as  DSMC  is  re¬ 
quired  simply  due  to  limited  resources.  It  is  currently  not  possible  to  calculate  the 
dynamics  of  each  molecule  for  a  large  system  using,  for  example,  Molecular  Dynamics 
(MD).  The  length  scales  in  MD  studies  involve  only  hundreds  of  molecules.  On  the 
other  hand,  the  number  of  molecules  on  the  anode  surface  of  a  magnetron  are  on  the 
order  of  10 15  per  square  centimeter.  Obviously,  MD  simulations  would  not  be  able 
to  sustain  such  a  large  number  of  molecules.  Hence  a  statistical  method  (DSMC)  is 
sought  which  makes  such  simulations  tractable. 

In  addition,  DSMC  techniques  are  being  considered  by  AFRL  for  HPM  simulations 
due  to  differences  in  time  scales  during  magnetron  operation  [156].  For  example, 
it  is  reasonable  for  a  magnetron  to  have  a  500-ns  pulse  width  with  0.1  s  between 
pulses.  The  duty  cycle  is  calculated  as  the  ratio  of  the  pulse  width  to  the  period.  In 
this  example,  the  duty  cycle  is  0.0005%.  In  other  words,  only  0.0005%  of  the  dwell 
time  is  occupied  with  a  pulse.  During  the  pulse  width,  high-energy  electrons  are 
explosively  emitted  from  the  plasma  created  on  the  cathode  surface.  Accurate  PIC 
methods  are  employed  at  the  pulse  time  scale  to  adequately  resolve  particle  positions 
and  velocities,  as  well  as  the  dynamic  electromagnetic  (EM)  field,  while  low-energy 
neutrals  can  be  assumed  to  be  stationary.  The  total  EM  field  is  comprised  of  an 
applied  portion  from  a  permanent  magnet,  and  a  dynamic  portion  created  by  the 
motion  of  the  electrons.  Once  high-velocity  particles  have  transferred  their  energy 
and  only  lower- velocity  particles  remain  in  the  cavity  (i.e.  between  pulses),  the  EM 
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field  may  be  considered  static,  and  DSMC  techniques  would  then  be  available  to 
describe  particle  interactions  on  a  larger  time  scale.  This  research  directly  supports 
their  efforts  by  providing  thermal  desorption  and  adsorption  models  for  the  anode 
surface  in  a  DSMC  framework. 

Adsorption  is  the  process  by  which  a  foreign  molecule  is  bound  to  a  surface. 
Absorption,  on  the  other  hand,  is  the  process  by  which  a  foreign  molecule  is  bound 
within  the  bulk.  Desorption  is  the  inverse  process  of  adsorption.  Therefore,  desorption 
is  the  process  by  which  the  bond  between  a  foreign  molecule  and  the  surface  is  broken. 

Some  of  the  most  elementary  surface  processes  in  important  plasma  systems  em¬ 
ployed  in  basic  research  and  in  modern  material  science  are  (dissociative)  adsorption 
and  (recombinative)  desorption  [27].  Modeling  these  processes  efficiently  will  allow 
for  more  accurate  predictions  of  system  performance  as  well  as  improved  system  de¬ 
signs  [7]. 

Unfortunately,  thermal  desorption  (and  adsorption)  dynamics  are  difficult  to  model. 
The  description,  interpretation,  and  prediction  of  state  and  energy  distributions 
of  desorbing  molecules  currently  represent  a  major  challenge  in  theoretical  chem¬ 
istry  [139].  It  should  come  as  no  surprise  that  the  understanding  of  surface  phe¬ 
nomena  is  still  far  from  that  of  gas-phase  reactions  [107],  since  a  gas-phase  reaction 
involves  only  a  small  number  of  molecules  when  compared  with  even  the  simplest  of 
gas-surface  systems. 

To  complicate  matters  further,  there  are  multiple  desorption  pathways  available, 
of  which  the  following  are  a  brief  sampling.  There  is  the  Langmuir-Hinshelwood 
(LH)  reaction,  characterized  by  two  adsorbed  atoms  combining  into  a  molecule  on 
the  surface  that  is  then  immediately  desorbed  [60;  146;  168],  and  the  Eley-Rideal 
(ER)  mechanism,  wherein  an  adatom  and  a  gas-phase  atom  combine  into  a  molecule 
that  is  desorbed  immediately  [28;  71;  133;  137].  In  the  presence  of  strong  EM  fields, 
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field-induced  desorption  (FID)  may  cause  adsorbed  molecules  to  desorb.  The  bulk 
material  itself  may  contribute  through  phonon- induced  desorption,  in  which  surface 
vibrations  sufficiently  excite  adatoms  to  desorb  [15].  In  addition,  collisions  from  ions, 
photons,  and  electrons  are  known  to  cause  desorption,  through  ion- induced  desorption 
(IID)  [21],  photo-induced  desorption  (PID)  [63],  and  ESD3  [4;  6;  41;  69;  70;  113;  121; 
140],  respectively. 

Some  of  these  desorption  pathways  are  present  in  a  typical  magnetron  system.  In 
the  magnetron  cavity,  a  plasma  forms  out  of  electrons,  ions,  and  neutral  molecules. 
Those  electrons  impact  the  anode  surface  after  being  emitted  from  the  cathode  and 
travelling  through  the  annulus.  Such  impacts  cause  either  ESD,  or  localized  surface 
heating  which  in  turn  induces  thermal  desorption  (LH).  When  the  ions  collide  with 
the  anode  surface,  they  either  scatter  or  induce  the  desorption  of  electrons,  ions, 
or  molecules  (HD).  Surface  coverage  may  repopulate  through  adsorption  from  ionic 
and  molecular  impacts,  as  well  as  from  permeation  and  diffusion  of  atoms  through 
the  bulk.  All  of  these  effects  combine  to  create  a  complex  andproblem  which  is 
currently  unresolved.  Questions  remain  about  the  relative  contributions  of  each  of 
these  effects,  questions  which  can  only  be  answered  with  experimental  and  theoretical 
contributions.  Such  work  will  allow  for  these  effects  to  be  understood  and  decoupled. 

Unfortunately,  all  of  the  relevant  surface  effects  cannot  be  addressed  in  this  re¬ 
search.  Therefore,  their  relative  contributions  must  be  prioritized.  The  two  most 
important  desorption  mechanisms  are  thought  to  be  caused  by  the  high  number  of 
electron  impacts  [59].  They  are  (1)  thermal  desorption  via  the  LH  mechanism,  in¬ 
directly  caused  by  localized  surface  heating  from  electron  impacts,  and  (2)  ESD. 
The  next  contributor  to  desorption  is  expected  to  be  thermal  desorption  via  the  ER 
mechanism,  indirectly  caused  by  localized  surface  heating  from  electron  impacts,  and 

3ESD  is  sometimes  referred  to  as  Desorption  Induced  by  Electronic  Transitions  (DIET). 
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directly  caused  by  ionic/molecular  impacts.  Finally,  adsorption  is  thought  to  occur 
from  ionic  and  molecular  surface  impacts. 

The  phenomena  of  desorption  and  adsorption  are  strongly  influenced  by  the  type 
of  chemical  bond  that  exists  between  the  adsorbate/gas  and  the  surface.  There  are 
two  classifications  available  in  the  literature,  physical  sorption  ( physisorption )  and 
chemical  sorption  ( chemisorption )  [62:178].  Physisorption  is  characterized  by  a  weak 
bond  dominated  by  van  der  Waals  forces,  usually  on  the  order  of  0.01  to  0.1  eV, 
and  does  not  involve  an  activation  barrier.  On  the  other  hand,  chemisorption  is 
characterized  by  a  much  stronger  bond,  typically  in  the  range  of  0.1  to  10  eV,  with 
the  inclusion  of  an  activation  barrier.  Chemisorption  is  the  result  of  either  hydrogen 
bonding,  covalent  chemical  bonding,  or  metallic  bonding. 

In  order  to  study  thermal  desorption  via  the  LH  pathway,  one  of  the  most  im¬ 
portant  surface  effects  in  HPM  magnetrons,  both  physisorptive  and  chemisorptive 
systems  will  be  investigated.  First,  the  inverse  process  of  desorption  (adsorption)  will 
be  considered  for  a  physisorption  system  in  Chapter  II.  This  system  will  be  as  simple 
as  possible,  where  the  adsorbate  is  xenon,  a  noble  gas,  and  the  surface  is  platinum. 
The  Xe-Pt  gas-surface  system  is  a  typical  starting  point  for  physisorption  modeling, 
and  will  provide  insight  into  how  to  model  surface  effects  in  a  DSMC  framework. 
Then  in  Chapter  III,  the  complex  process  of  thermal  desorption  (LH)  will  be  studied 
on  the  simplest  chemisorption  gas-surface  system  available,  H2-Cu  [126].  Fortunately, 
H2-C11  is  also  the  most  thoroughly  studied  example  of  activated  adsorption  [114],  The 
study  of  H2-Cu  desorption  will  provide  both  simulation  results  directly  applicable  to 
current  HPM  research,  and  a  new  methodology  from  which  desorption  can  be  modeled 
for  any  gas-surface  system  in  the  future. 
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1.1  Physisorption 


Physisorption  is  investigated  (adsorption  and  scattering)  in  Chapter  II  for  a  model 
system,  Xe-Pt(lll).  This  gas-surface  combination  is  typically  employed  by  researchers 
when  developing  physisorption  models.  Xenon  is  a  noble  gas,  therefore  it  adsorbs  and 
desorbs  monatomically.  There  are  no  internal  energies  to  consider.  There  is  no  ac¬ 
tivation  energy  and  the  adsorbate  bond  is  weak.  Also,  the  Potential  Energy  Surface 
(PES)  is  relatively  smooth. 

Even  though  physisorption  is  modeled  only  for  adsorption  and  scattering,  all  three 
surface  processes  (adsorption,  scattering,  and  desorption)  are  related.  Adsorption 
and  scattering  have  in  common  the  collision  of  a  gas  molecule  with  the  surface.  The 
link  between  adsorption  and  desorption  is  even  more  direct,  because  they  are  inverse 
processes  of  one  another  [78:161-166].  The  former  is  concerned  with  a  gas-phase 
molecule  adhering  to  the  surface,  while  the  latter  considers  a  molecule  on  the  surface 
being  released  to  the  gas-phase.  An  understanding  of  one  of  these  processes  invites 
greater  insight  into  the  other. 

The  assumption  that  the  problem  is  symmetric  with  respect  to  time  is  termed  mi¬ 
croscopic  reversibility,  a  stricter  form  of  the  principle  of  detailed  balance.  Time  sym¬ 
metry,  or  T-symmetry,  indicates  that  the  physical  laws  are  symmetric  whether  time 
is  advanced  forwards  or  backwards.  Detailed  balance  is  only  applicable  under  equilib¬ 
rium  conditions.  However,  microscopic  reversibility  applies  even  in  non-equilibrium 
situations.  The  principle  of  detailed  balance  provides  a  mechanism  by  which  the  prop¬ 
erties  of  a  system  in  equilibrium  are  maintained  [29].  Population  and  depopulation 
rates  of  every  microscopic  and  macroscopic  state  are  equal  under  detailed  balance. 
However,  microscopic  reversibility  assumes  that  the  trajectory  and  state  of  a  molecule 
can  be  simulated  forwards  and  backwards  in  time  with  identical  results,  which  holds 
under  Hamiltonian  dynamics.  Thus,  under  microscopic  reversibility,  a  molecule  des- 
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orbing  from  a  surface  provides  information  on  a  molecule  adsorbing  to  a  surface  under 
identical  conditions.  Recently,  microscopic  reversibility  was  shown  to  also  hold  for 
open  systems,  bridging  the  first  and  second  laws  of  thermodynamics  [119]. 

The  novel  adsorption  and  scattering  DSMC  method  developed,  called  the  Mod¬ 
ified  Kisliuk  with  Scattering  (MKS)  method,  is  simple.  The  main  points  are:  (1) 
determine  if  an  impinging  gas  molecule  collides  with  an  adatom  or  with  the  surface 
by  comparing  a  random  uniform  number  with  the  fractional  surface  coverage,  (2) 
calculate  the  molecule’s  scattering  energy  and  angles  by  incorporating  an  existing 
scattering  computational  method  [Cercignani-Lampis-Lord  (CLL)],  (3)  compare  the 
scattering  energy  with  the  well-depth  parameter  of  the  PES  to  determine  whether 
the  molecule  adsorbs  or  scatters,  and  (4)  repeat  as  necessary. 

In  spite  of  its  simplicity,  MKS  recreates  experimental  data  for  Xe-Pt(lll)  remark¬ 
ably  well.  Thus,  a  strength  of  the  DSMC  method  is  that  underlying  assumptions  can 
be  included  in  the  model  in  a  relatively  simple  manner.  Analytical  expressions  do 
not  need  to  be  developed  every  time  an  assumption  is  added  or  altered.  One  can  pick 
and  choose  the  assumptions  deemed  relevant,  and  build  on  them  as  more  complex 
problems  are  considered.  MKS  is  actually  a  transformation  (or  translation)  of  an 
existing  model  [Modified  Kisliuk  (MK)]  into  a  DSMC  algorithm,  although  this  was 
not  recognized  originally  as  MKS  was  being  developed.  MKS  was  built  independently 
from  MK  with  assumptions  that  expressed  the  dominant  underlying  physical  phenom¬ 
ena  of  the  Xe-Pt(lll)  system.  Once  MKS  results  were  matched  with  experimental 
data,  it  was  discovered  that  the  assumptions  behind  MKS  perfectly  coincided  with 
the  assumptions  that  went  into  MK.  Hindsight  thus  showed  that  analytical  models 
that  cannot  be  directly  utilized  in  DSMC  applications,  can  however  be  translated 
into  a  DSMC  algorithm  simply  by  implementing  their  underlying  assumptions.  Fur¬ 
thermore,  existing  models  can  be  augmented,  or  extended,  by  this  method.  The  MK 


model  only  calculates  the  probability  that  a  gas  molecule  adsorbs  as  a  function  of 
current  adsorbate  coverage.  However,  in  addition  to  the  adsorption  probability,  MKS 
also  calculates  the  scattering  energy  and  angles. 

1.2  Chemisorption 

Thermal  desorption  (LH)  is  modeled  in  Chapter  III  for  H2  on  the  (100)  surface 
of  copper  [H2-Cu(100)],  a  chemisorptive  system.1  Hydrogen  is  the  most  prevalent 
residual  gas  in  metal  vacuum  systems  at  very  low  pressures  [Ultra-High  Vacuum 
(UHV)  and  Extreme-High  Vacuum  (XHV)].  In  fact,  researchers  who  seek  to  lower  the 
pressure  of  a  metal  system  down  to  10” 10  Pa  (XHV)  view  the  mitigation  of  hydrogen 
outgassing  as  the  most  challenging  problem  they  face  [141],  and  is  investigated  in  the 
literature  [40;  42;  43;  171].  The  H2-Cu  gas-surface  system  is  considered  because  it  is 
the  most  commonly-studied  example  of  activated  chemisorption.  The  (100)  face  of 
copper  was  chosen  simply  because  there  is  sufficient  modeling  and  experimental  data 
in  the  literature  for  the  H2-Cu(100)  system.  Some  other  resources  in  the  literature 
include  two  methods  of  doping  hydrogen  onto  copper  for  desorption  studies  [5;  51], 
and  various  investigations  into  thermal  desorption  (LH)  for  the  D2-Cu  [115;  116]  and 
H2-Cu  [25;  56;  87;  132;  144;  159]  systems. 

Hydrogen  atoms  tend  to  permeate  metals  relatively  easily  [54;  55;  130].  Thus,  even 
if  a  surface  is  completely  clear  of  adsorbate,  hydrogen  desorption  and  outgassing  may 
still  occur  during  system  operation.  Sometimes,  the  hydrogen  will  first  permeate, 
then  bond  in  sub-surface  sites  in  the  metal  (even  though  there  is  some  discussion 
on  this  matter)  [102;  124].  For  hydrogen  permeating  copper,  it  is  possible  that  H2 
recombines  in  a  subsurface  state  [114],  and  receives  additional  energy  as  it  desorbs 
due  to  the  (bulk)  absorption  energy  barrier  [33] .  Typically  in  desorption  experiments, 
4The  nomenclature  for  surface  structure  [i.e.  (100)]  is  discussed  in  Section  3.1. 
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hydrogen  is  supplied  to  the  surface  via  thin-membrane  permeation.  However,  one  can 
supply  hydrogen  directly  to  the  surface  via  an  atomic  beam  in  order  to  avoid  any 
complicating  permeation  effects  [146]. 

LH  reactions  are  driven  by  the  thermal  energy  of  the  surface  being  transferred  to 
the  adsorbate;  they  are  rarely  caused  by  adsorbate  excitations,  since  the  energy  trans¬ 
fer  rate  from  the  lattice  to  the  adsorbate  is  typically  13  orders  of  magnitude  greater 
than  the  reaction  rate  [60].  However,  LH  desorption  events  result  in  non-equilibrium 
distributions,  heavily  influenced  by  the  curvature  of  the  PES  [168].  One  reason  for 
this  result  is  that  the  potential  energy  barrier  to  recombination  at  the  surface  accel¬ 
erates  desorbing  molecules  away  from  the  bulk  (and  inhibits  adsorption)  [146],  thus 
providing  desorbed  products  with  higher-than-average  kinetic  energies  than  would 
otherwise  be  available  simply  from  thermal  equilibrium  with  the  surface.  In  addi¬ 
tion,  the  desorption  of  hydrogen  from  copper  is  an  early-barrier  process,  implying 
that  strong  energy  transfers  occur  after  recombination  [139],  further  perturbing  the 
molecule’s  distributions  from  equilibrium. 

Even  though  it  will  not  be  studied  here,  activated  dissociative  chemisorption  (ad¬ 
sorption)  is  an  important  area  of  study  since  it  is  the  rate-limiting  step  in  many 
gas-surface  reactions  [145].  It  is  mostly  a  function  of  the  normal  component  of  the  in¬ 
cident  energy,  a  relationship  known  as  normal  energy  scaling  [5] .  However,  the  surface 
temperature  also  aids  in  molecular  dissociation  [120].  Also,  adsorption  is  enhanced 
for  vibrationally  excited  H2  molecules,  where  around  60%  of  the  vibrational  energy 
is  expended  in  overcoming  the  activation  barrier  [114],  In  order  to  understand  how 
each  of  these  mechanisms  contributes  to  adsorption,  state-resolved  adsorption  [76]  and 
scattering  [169;  170]  has  been  studied  for  H2-Cu(100).  Electronically  non-adiabatic 
effects  are  thought  to  affect  scattering,  where  hydrogen  molecules  lose  some  vibra¬ 
tional  energy  while  exciting  copper  electrons  [101].  Other  areas  of  interest  include 
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low-temperature  adsorption  [H2-Cu(100)]  [138],  how  to  prepare  a  copper  surface  for 
scattering  studies  [12],  and  various  theoretical  and  experimental  publications  [H2- 
Cu]  [8;  34;  36;  38;  56;  58;  85;  90;  123;  158;  159;  161;  162], 

Unfortunately,  chemisorption  cannot  successfully  be  modeled  in  the  same  man¬ 
ner  as  physisorption.  Physisorption  concerns  itself  with  systems  in  which  weak 
van  der  Waals  interactions  dominate,  ffowever,  in  chemisorptive  systems,  the  elec¬ 
tronic  structure  of  the  adsorbate  is  significantly  modified  due  to  the  presence  of  the 
surface.  As  such,  the  corrugation  and  shape  of  the  PES  can  play  a  large  role  in 
chemisorptive  dynamics.  Thus,  a  different  approach  is  followed  for  chemisorption 
modeling,  which  is  to  model  the  trajectory  of  each  desorbing  molecule  in  a  Classical 
Trajectory  (CT)  formulation  under  the  influence  of  a  six-dimensional  (6D)  PES. 

The  PES  describes  the  interaction  between  a  surface  and  a  molecule.  From  the 
PES,  one  can  calculate  the  forces  acting  on  the  molecule,  and  given  appropriate  initial 
conditions,  the  molecule’s  position  and  velocity  as  well.  For  a  single  atom,  only  three 
degrees  of  freedom,  or  three  dimensions  (3D),  are  required  to  adequately  model  the 
PES.  For  a  diatomic  molecule,  six  degrees  of  freedom  are  needed.  Likewise,  for  a 
molecule  with  N  atoms,  there  are  potentially  3N  degrees  of  freedom  required  [65]. 
One  may  simplify  the  calculations  by  neglecting  some  of  the  degrees  of  freedom, 
ffowever,  it  has  been  shown  that  all  six  dimensions  must  be  considered  when  modeling 
associative  (dissociative)  chemisorption  for  a  diatomic  molecule  [44;  45].  Even  though 
all  six  degrees  of  freedom  are  addressed  here,  there  is  assumed  to  be  no  change  in 
the  electronic  state  of  the  system  during  the  reaction.  In  other  words,  the  reactions 
are  assumed  to  be  electronically  adiabatic.  Also,  since  the  motion  of  the  admolecule 
can  be  obtained  from  the  PES,  it  is  not  necessary  to  identify  frustrated  (or  hindered) 
degrees  of  freedom.5 

5Since  adsorbed  molecules  behave  differently  than  free  molecules,  their  motions  are  sometimes 
referred  to  as  frustrated  or  hindered,  such  as  hindered  rotation  or  hindered  translation  [17;  26;  66; 


II 


Lattice  motion  is  assumed  to  be  insignificant.  Thus,  a  Rigid  Surface  (RS)  as¬ 
sumption  is  employed,  which  considers  bulk  atoms  frozen  in  time  and  space  [60].  One 
reason  for  this  assumption  is  that  phonons  can  safely  be  neglected  when  modeling 
surface  effects  for  a  light  molecule  such  as  H2  [1;  2;  122;  149;  150;  151;  158].  Also, 
the  time  scale  for  hydrogen  association  (dissociation)  on  copper  is  much  smaller  than 
that  of  the  substrate  motion.  Another  justification  is  that  surface  reconstruction  is 
unimportant  in  this  case,  since  adsorbed  hydrogen  perturbs  the  substrate  lattice  only 
minimally  [60].  For  a  more  refined  look,  however,  the  reader  is  referred  to  current 
investions  on  how  small  perturbations  in  lattice  position  away  from  equilibrium  can 
effect  the  potential  energy  landscape  of  the  H2-Cu  interaction  [22;  177]. 

Wiesenekker  et  al.  [83;  84;  174]  developed  a  6D  PES  to  describe  the  molecule- 
surface  interaction  in  the  H2-Cu(100)  system.  This  PES  is  based  on  slab  calculations 
with  the  Generalized  Gradient  Approximation  (GGA)  in  DFT.  It  is  expressed  as  a 
fit  to  an  analytical  form,  bypassing  the  need  for  computationally  intensive  on-the- 
fly  DFT  calculations.  An  analytical  fit  also  provides  a  high  degree  of  flexibility  in 
matching  the  PES  to  experimental  and  simulated  energy  values.  The  full  6D  PES  is 
achieved  by  splicing  together  eight  two-dimensional  (2D)  PES  slices  [175;  176],  which 
are  individually  anchored  to  experimental  results.  This  6D  PES  was  specifically 
tailored  to  be  accurate  in  the  region  between  the  gas  phase  and  the  Transition  State 
(TS)  (the  location  or  zone  where  the  individual  adatoms  combine  into  a  molecule). 
Within  the  TS,  the  PES  is  fairly  accurate,  and  between  the  TS  and  the  surface,  the 
PES  is  inaccurate. 

There  are  other  options  when  considering  a  PES  for  the  H2-Cu(100)  system  [32;  37; 
48].  One  is  the  common  London-Eyring-Polyani-Sato  (LEPS)  form  [109;  154],  known 
for  its  computational  simplicity  and  speed.  The  LEPS  form  assumes  that  atom-atom 
110;  129;  178], 
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and  atom-surface  interactions  are  adequately  described  by  Morse  potentials  [109]. 
However,  the  LEPS  form  contains  only  a  few  adjustable  parameters,  and  is  therefore 
less  flexible  than  the  Wiesenekker  PES.  Another  option  is  the  PES  constructed  by 
Olsen  et  al.  [126].  Unfortunately,  they  did  not  publish  their  fitting  parameters,  and 
hence  their  PES  cannot  be  reproduced.  It  should  be  noted  that  the  methods  presented 
here  are  applicable  to  any  PES  (LEPS,  analytical,  etc.).  In  fact,  studies  concerning 
H2-Cu  PESs  continue  to  be  reported. 

The  trajectories  of  desorbing  molecules  are  modeled  within  the  CT  formulation. 
CT  simulations  ignore  the  quantum  effects  of  tunneling  and  interferences  [88;  118]. 
These  effects  may  be  neglected  for  relatively  heavy  molecules,  and  for  surfaces  at 
temperatures  above  their  Debye  temperature  [167].  The  Debye  temperature  of  a 
solid  is  interpreted  as  the  temperature  at  which  all  phonon  modes  are  excited.  Even 
though  H2  is  considered  to  be  a  light  molecule,  the  temperatures  considered  in  this 
work  for  copper  are  well  above  its  Debye  temperature  of  344  K.  Therefore,  the  CT 
formulation  is  appropriate  here.  A  significant  benefit  of  modeling  desorption  with 
CT  simulations  is  that  rovibrational  state  distributions  can  be  calculated,  which  have 
been  shown  to  be  important  in  H2-Cu  dynamics  [152],  As  an  aside,  the  Quasi-Classical 
Trajectory  (QCT)  method  attempts  to  include  quantum  effects  by  incorporating  Zero- 
Point  Energy  (ZPE)  in  the  initial  conditions  [19;  108]. 

Desorption  calculations,  in  order  to  be  incorporated  into  DSMC  codes,  need  to 
provide  information  about  the  molecule  including  its  velocity,  and  its  rotational  and 
vibrational  energies.  Fortunately,  CT  simulations  are  able  to  do  just  that.  Recently, 
kinetic  Monte  Carlo  (MC)  simulations  of  surface  reactions  have  been  coupled  with 
continuum  computational  fluid  dynamics  (CFD)  methods  [155].  However,  there  is  no 
indication  in  the  literature  that  CT  has  been  developed  for  use  in  DSMC  simulations 
before  this  work. 
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Also,  for  time-accurate  simulations,  a  method  for  timing  desorption  events  must 
be  included.  Timing  will  be  incorporated  by  noting  that  each  desorption  event  for 
H2-Cu(100)  can  be  assumed  to  be  independent  of  one  another,  thus  representing  a 
Poisson  process.  A  rate  of  surface  activity  then  allows  for  the  time  between  successive 
events  to  be  calculated.  In  this  manner,  desorption  events  will  occur  during  specific 
times  in  the  DSMC  simulation.  More  about  this  topic  will  be  discussed  in  Section  3.8. 

1.3  Contributions 

In  summary,  this  research  develops  thermal  desorption  and  adsorption  models  in 
a  framework  to  incorporate  into  Direct  Simulation  Monte  Carlo  codes.  As  a  result, 
plasma-surface  modeling  will  greatly  benefit.  Current  state-of-the-art  plasma-surface 
interaction  modeling  does  not  capture  essential  molecular  physics.  For  example,  ad¬ 
sorption  is  represented  by  a  single  probability,  scattering  is  calculated  with  a  proba¬ 
bility  curve,  and  desorption  is  assumed  to  exhibit  equilibrium  distributions  [165]. 

This  work  is  a  significant  step  forward  for  the  aeronautical  and  mechanical  en¬ 
gineering  communities  in  modeling  dynamic  surface  phenomena  (both  physisorp- 
tive  and  chemisorptive)  at  the  microscopic  level.  These  fields  have  not  previously 
had  access  to  microscopic  surface  modeling  within  the  Direct  Simulation  Monte 
Carlo  framework.  Microscopic  methods  have  been  neglected  in  favor  of  macro¬ 
scopic  fluid  approaches  focusing  on  reaction  rates  [103;  105;  163].  Unfortunately, 
macroscopic  approaches  are  not  able  to  model  in  detail  the  state-specific  dynamics 
of  desorbing  molecules,  even  when  considering  many  adsorbate  spatial  configura¬ 
tions  [47;  50;  67;  68;  77;  79;  80;  81;  82;  91;  100;  104;  112;  125;  128;  173;  181;  183], 

In  Chapter  II  (Physisorption),  a  new  algorithm  is  developed,  called  the  Modi¬ 
fied  Kisliuk  with  Scattering  method.  This  model  is  the  first  of  its  kind  to  account 
for  both  adsorption  and  scattering.  It  accurately  calculates  adsorption  probabili- 
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ties.  It  also  predicts  super-elastic  scattering  profiles,  a  capability  currently  lacking 
in  Direct  Simulation  Monte  Carlo  models.  Scattering  calculations  are  even  improved 
from  state-of-the-art  models.  Gas-adlayer  interactions  are  considered,  a  first  for  the 
Cercignani-Lampis-Lord  kernel.  Finally,  it  is  demonstrated  how  the  Modified  Kisliuk 
with  Scattering  method  can  be  used  to  determine  accommodation  coefficients. 

In  Chapter  III  (Chemisorption),  another  novel  method  is  introduced  to  model 
thermal  desorption  for  Direct  Simulation  Monte  Carlo  applications.  The  equations  of 
motion  are  non-dimensionalized,  a  description  surprisingly  not  found  in  the  literature. 
Multiple  transition  states  are  incorporated  into  the  calculations,  whereas  researchers 
in  the  held  have  only  considered  a  single  transition  state  until  now.  A  truncated 
Maxwell- Boltzmann  distribution  is  developed  and  its  accept-reject  form  is  derived. 
Correct  event  timing,  adapted  here  for  Direct  Simulation  Monte  Carlo  applications, 
is  included  for  temporally-accurate  simulations.  It  is  shown  that  Gaussian  weighting 
is  not  required  for  this  method,  an  unexpected  result  which  greatly  reduces  compu¬ 
tational  run  times.  Also,  the  absorption  barrier  energy  due  to  permeation  and/or 
diffusion  is  shown  to  significantly  contribute  to  thermal  desorption.  Specifically,  it 
directly  affects  the  translational  energy,  while  it  has  little  impact  on  the  internal 
energies  of  desorbing  molecules. 
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II.  Physisorption 


There  is  an  immediate  opportunity  to  improve  surface  effects  modeling  in  the 
DSMC  field.  Currently,  adsorption  is  completely  ignored  due  to  its  complex  nature. 
The  modeling  of  scattering  received  significant  attention  a  little  less  than  20  years  ago. 
However,  there  has  been  little  improvement  since  then.  The  simplest  assumption  for 
scattering  is  that  molecules  do  so  specularly.  HPM  simulations  usually  only  include 
this  specular  reflection  boundary  condition.  Other  fields,  such  as  hypersonics,  have 
adopted  a  more  advanced  technique  with  the  CLL  kernel;  the  CLL  boundary  condition 
is  capable  of  producing  a  somewhat  realistic  lobular  scattering  profile  [97]. 

This  chapter  details  the  development  and  results  of  a  new  model,  MKS,  that  ad¬ 
dresses  both  adsorption  and  scattering  for  a  physisorptive  system  in  a  DSMC  frame¬ 
work.  First,  a  brief  introduction  to  simple  adsorption  models  will  be  given,  building 
to  a  discussion  of  the  MK  model.  Then,  the  MKS  method  will  be  developed,  including 
a  detailed  discussion  of  the  CLL  kernel.  Finally,  adsorption  results  from  MKS  will 
be  compared  with  experimental  data  with  respect  to  the  Xe-Pt(lll)  gas-surface  sys¬ 
tem,  illustrating  the  benefits  of  MKS.  Besides  being  the  first  method  to  handle  both 
adsorption  and  scattering  in  a  DSMC  framework,  MKS  also  improves  the  scattering 
properties  of  the  state-of-the-art  DSMC  scattering  model. 

2.1  Previous  Adsorption  Models 

This  section  introduces  simple  adsorption  models,  in  preparation  to  discuss  the 
new  MKS  model  in  Section  2.2.  Starting  with  Henry’s  Law,  the  historical  models  will 
progressively  increase  in  complexity,  until  an  adsorption  model  capable  of  handling 
the  behavior  of  Xe-Pt(lll)  is  explained  (MK). 
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2.1.1  Henry’s  Law. 


The  simplest  adsorption  model,  Henry’s  Law,  assumes  (1)  a  gas  molecule  will  ad¬ 
sorb  independently  of  the  adlayer,  and  (2)  all  adsorption  sites  are  equivalent  [62:175- 
197].  For  this  model,  the  adsorption  probability  is  independent  of  coverage,  or 


S(9)  =  S0, 


(1) 


where  6  is  the  adlayer  coverage  and  So  is  the  initial  adsorption  probability.1 

2.1.2  Langmuir  Model. 

The  Langmuir  model  is  slightly  more  complicated  than  Henry’s  Law  because  it 
assumes  (1)  adatoms  occupy  specific  adsorption  sites  on  the  surface,  (2)  there  is  no 
adatom  mobility,  (3)  a  gas  molecule  does  not  adsorb  if  it  encounters  an  occupied 
site,  and  (4)  there  are  no  lateral  interactions  [62:240;  106].  This  model  predicts  that 
a  gas  molecule  is  adsorbed  into  an  unoccupied  site  with  probability  S0,  while  gas 
molecules  are  never  adsorbed  onto  occupied  sites.  Thus,  the  adsorption  probability 
S(9)  decreases  linearly  with  coverage, 


S(6)  =  S0(l-6). 


(2) 


2.1.3  Kisliuk  Model. 

Kisliuk  developed  a  model  that  considered  adsorption  via  both  an  intrinsic  and  an 
extrinsic  precursor,  with  no  direct  Langmuirian  adsorption  [10].  An  intrinsic  precursor 
in  adsorption  is  a  molecule  that  is  occupying  an  adsorption  site,  but  that  has  not  fully 
adsorbed  to  the  surface.  In  other  words,  the  molecule  is  still  mobile  on  the  surface. 

■'Tor  a  full  listing  of  all  the  symbols  used  in  this  document,  including  their  respective  units,  please 
refer  to  Appendix  A. 
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Similarly,  an  extrinsic  precursor  is  a  molecule  that  occupies  an  adsorption  site  above 
the  adlayer,  which  has  not  fully  adsorbed  to  the  adlayer  (or  surface).  In  other  words, 
precursors  are  transitory  states  of  an  adatom  as  it  either  adsorbs  or  desorbs.  In 
steady- state, 


S{0) 


St) 


1  +  K 


e  ’ 


1-9 


(3) 


where  So  and  K  are  functions  of  the  trapping  probability  into  precursor  states  as  well 
as  rate  constants  for  the  intrinsic  and  extrinsic  precursor  states.  This  model  assumes 
(1)  no  lateral  interactions  and  (2)  a  random  adlayer  configuration.  It  turns  out  that 
according  to  this  model,  S(6)  also  decreases  with  increasing  coverage.  However,  there 
are  systems  where  the  opposite  trend  holds,  such  as  Xe-Pt(lll)  [10]. 


2.1.4  Modified  Langmuir  Model. 

Another  model,  which  will  here  be  referred  to  as  a  modified  Langmuir  model, 
was  proposed  to  fit  experimental  data,  which  assumed  that  the  adsorption  pathways 
included  both  direct  Langmuirian  adsorption  and  mobile  adsorption  onto  the  second 
adlayer.  In  other  words, 

S(9)  =  S0(l  -9)  +  S*9 ,  (4) 

where  Sq  is  the  adsorption  probability  onto  the  adlayer.  For  Sq  >  So,  S(9)  increases 
with  coverage. 


2.1.5  Modified  Kisliuk  Model  (MK). 

The  MK  model  was  developed  by  combining  the  original  Kisliuk  model  with  the 
modified  Langmuir  model  [10].  MK  assumes  (1)  Langmuirian  adsorption  for  a  gas 
molecule  in  contact  with  an  unoccupied  site  with  adsorption  probability  So,  (2)  ex¬ 
trinsic  precursor  mobile  adsorption  with  adsorption  probability  Sq,  (3)  the  configu- 
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ration  of  the  adlayer  is  random,  and  (4)  no  lateral  interactions  exist.  In  steady-state, 
assuming  that  the  extrinsic  precursor  coverage  is  small, 


S(0)  =  S0(l  -  9)  +  S0*(l  -  (5) 

t  QmU 

where  qm  is  a  function  of  extrinsic  precursor  rate  constants. 

2.2  Novel  Adsorption  and  Scattering  Method 

With  that  brief  background  in  mind,  a  novel  adsorption  and  scattering  method 
is  developed  in  this  section.  Termed  the  Modified  Kisliuk  with  Scattering  (MKS) 
method,  this  new  technique  builds  on  the  assumptions  of  MK  for  DSMC  and  the 
Multi-Stage  (MS)  gas-surface  interaction  model.  It  is  capable  of  modeling  both  ad¬ 
sorption  and  scattering,  an  accomplishment  previously  unattained.  Further,  scatter¬ 
ing  properties  of  MKS  are  shown  to  be  superior  to  those  of  state-of-the-art  calcu¬ 
lations  with  the  CLL  kernel.  MKS  assumes  (1)  no  precursor- mediated  adsorption 
for  impingement  on  an  unoccupied  site,  (2)  mobile  precursor-mediated  adsorption 
for  impingement  on  the  adlayer,  (3)  the  adlayer  configuration  is  random,  and  (4)  no 
lateral  interactions  are  present. 

2.2.1  Multi-Stage  (MS)  Gas-Surface  Interaction  Model. 

The  implementation  of  MKS  is  similar  to  that  of  the  MS  gas-surface  interaction 
model  [180].  The  MS  model  is  applicable  to  the  thermal  scattering  regime  [53:103], 
where  the  interactions  produce  lobular  and  diffuse  scattering.  The  MKS  method 
presented  in  this  work  is  also  limited  to  this  regime,  and  one  should  be  aware  of  the 
conditions  under  which  it  loses  applicability.  The  radius  parameter  1Z  is  the  main 
determining  factor  of  the  scattering  regime,  and  a  value  of  1Z  >  1.5  indicates  thermal 
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scattering.  1Z  is  defined  as 


(6) 


where  Rcg  is  the  effective  interaction  radius,  or  the  smallest  distance  between  the  gas 


and  surface  atoms  during  a  gas-surface  collision,  and  Rc  is  a  critical  value  of  RcS.  For 


a  given  incident  translational  energy  Etr,  Rcs  is  estimated  from  the  one-dimensional 
(ID)  approximation 


where  eLJ  is  the  Lennard-Jones  (LJ)  well-depth  parameter  for  the  interaction  between 
the  gas  molecule  and  a  surface  atom,  and  is  the  equilibrium  distance  between  the 
gas  molecule  and  a  surface  atom.  Rc  is  then  the  distance  from  the  gas  molecule  to  a 
surface  atom  as  it  crosses  the  surface  plane.  For  the  commensurate  layer  on  a  (111) 


surface,  the  critical  value  is  Rc  =  do/\/3,  where  do  is  the  equilibrium  distance  between 


surface  atoms.2 

The  MS  model  breaks  down  the  gas-surface  collision  into  three  steps.  First,  trans¬ 
lational  and  rotational  energies  are  calculated  from  a  model  equation,  which  is  based 
on  MD  simulations  of  that  specific  collision  with  those  specific  gas  and  surface  species. 
Next,  from  the  PES,  the  scattering  angles  are  found.  Finally,  depending  on  the  re¬ 
sults  from  the  first  two  steps,  the  gas  molecule  either  scatters,  re-collides  with,  or 
is  trapped  by  the  surface.  If  the  gas  molecule  re-collides  with  the  surface,  then  the 
process  is  repeated,  until  a  maximum  of  ten  collisions  occur.  At  that  point,  the  gas 
molecule  is  trapped  by  the  surface,  and  is  assumed  to  desorb  diffusely  at  the  surface 
temperature  Ts  during  that  same  time  step.  Trapping  is  defined  for  an  adatom  in 
a  weakly-bound  mobile  state.  However,  once  the  adatom  is  more  strongly  bound  to 
the  surface,  it  is  said  to  stick  or  adsorb.  The  MS  model  only  calculates  scattering 

2The  commensurate  layer  on  a  (111)  surface  is  also  described  as  (vdlx  ^/3)R  30°.  The  unit  surface 
cell  of  adsorbate  is  configured  in  a  rhombus,  with  all  four  sides  of  length  \/3  times  the  length  of  the 
bulk  unit  cell,  and  rotated  30°  from  the  bulk  unit  cell’s  orientation. 
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and  trapping  over  a  simulation  time  step,  meaning  that  molecules  trapped  on  the 
surface  are  scattered  by  the  end  of  that  time  step;  it  does  not  consider  adsorption.  In 
fact,  the  literature  does  not  appear  to  contain  any  DSMC  adsorption  (or  desorption) 
models. 

2.2.2  New  Algorithm. 

The  MKS  method,  on  the  other  hand,  can  handle  adsorption  as  well  as  scattering. 
Like  MS,  MKS  also  breaks  up  the  gas-surface  collision  into  stages: 

1.  Determine  if  the  gas  molecule  collides  with  an  adatom  or  with  the  surface.  A 
uniform  random  number  0  <  R  <  1  is  compared  with  the  adlayer  coverage 
0  <  0  <  1.  If  R  <  0,  then  the  gas  molecule  collides  with  an  adatom.  Otherwise, 
the  gas  molecule  collides  with  the  surface.  6  is  measured  in  monolayers  (ML). 

2.  The  gas  molecule  scatters  using  a  scattering  method  or  kernel.  Accommoda¬ 
tion  coefficients,  defined  below,  are  provided  for  both  the  gas-surface  ( an ,  at) 
and  the  gas-adatom  (anjad,  ott, ad)  systems,  ct^ad  and  cpiad  are  new  coefficients, 
introduced  with  MKS  to  account  for  the  gas-adlayer  interactions.  Note  that 
Langmuirian  kinetics  are  recovered  for  cp^d,  cp,ad  =  0. 

3.  If  the  post-collisional  normal  translational  energy  E[v  n  of  the  gas  molecule  is 
sufficient  to  escape  the  PES  {E[r  n  >  2eLJ),  then  it  scatters.  Here,  eLJ  is  the 
Lennard-Jones  potential  well  parameter  for  the  interaction  between  the  gas 
molecule  and  a  surface  molecule  ( not  the  entire  surface).  The  threshold  value 
of  2eLJ  was  determined  by  Yamanishi  et  al.  and  is  discussed  below  [180]. 

4.  If  the  gas  molecule  cannot  escape  the  PES,  but  has  enough  total  energy  E'tot  to 
overcome  the  adsorption  energy  and  has  not  collided  with  the  surface  more  than 
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ten  times,  then  it  re-collides  with  the  surface.  In  other  words,  if  E'tot  >  E^s 
and  Nco\\s  <  10,  then  return  to  Step  1  using  post-collisional  values. 

5.  Otherwise,  the  gas  molecule  is  adsorbed. 

The  MKS  method  differs  from  the  MS  model  in  a  few  respects.  First,  any  scatter¬ 
ing  method  or  kernel,  including  a  model  equation,  can  be  used  to  determine  scattering 
values  in  MKS,  while  MS  requires  a  tailor-made  model  equation  for  each  separate  case. 
The  CLL  scattering  kernel  [97]  was  chosen  since  current  DSMC  codes  already  use  the 
CLL  kernel.  However,  MKS  is  modular  in  that  any  scattering  kernel  or  model  could 
be  employed  in  the  place  of  the  CLL  kernel.  CLL  requires  accommodation  coeffi¬ 
cients  ( an ,  at,  aniad,  and  attSL d)  to  be  specified,  and  so  they  are  utilized  here.  But  if 
a  different  scattering  kernel  or  method  were  to  be  incorporated  into  MKS  instead  of 
CLL,  then  different  parameters  would  be  required.  A  contribution  made  here  is  ap¬ 
plying  the  CLL  kernel  to  adlayer  scattering.  Up  until  now,  only  surface  scattering  has 
been  considered  with  CLL.  Next,  the  MS  model  considers  only  a  clean  surface  while 
MKS  accounts  for  an  adlayer  being  present,  up  to  one  monolayer.  This  limitation 
illustrates  a  drawback  of  MS,  because  the  model  equation  must  be  developed  for  each 
individual  gas-surface  system.  In  fact,  adsorbates  alter  the  PES  of  the  interaction, 
and  therefore  a  multitude  of  model  equations  would  ideally  have  to  be  developed  for 
a  single  gas-surface  system.  The  new  MKS  model  is  not  limited  to  perfectly-clean 
surfaces.  Additionally,  the  criteria  have  been  changed  for  determining  whether  or 
not  a  gas  molecule  overcomes  the  potential  well  of  the  PES.  From  MD  simulations, 
Yamanishi  et  al.  claim  that  when  a  gas  molecule  scatters  off  of  a  surface,  it  has  an 
average  potential  energy  of  2eLJ  [180].  In  the  MS  model,  the  gas  molecule  scatters  if 
both  the  post-collisional  normal  and  tangential  translational  energies  are  greater  than 
or  equal  to  this  average  potential  energy  (i.e.,  E[r  n  >  2eLJ  and  E[r  f  >  2 eLJ).  However, 
while  comparing  MKS  with  experimental  adsorption  probability  data,  this  research 
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found  that  the  gas  molecule  should  scatter  only  if  E[r  n  >  2eLJ.  Hence,  no  restriction 
was  placed  on  E[rt  in  MKS.  A  sensitivity  analysis  was  performed  to  determine  the 
optimal  threshold  value.  It  was  found  in  this  work  that  2eLJ  is  reasonable.  A  brief 
discussion  of  this  finding  is  in  Section  2.4.  Finally,  MS  includes  rotational  energy 
while  MKS  in  its  current  form  does  not  consider  internal  energy  modes. 


2.2.3  Accommodation  Coefficients. 


Before  the  typical  CLL  kernel  can  be  reviewed,  accommodation  coefficients  should 
be  discussed  since  the  CLL  kernel  requires  them  as  inputs.  Accommodation  coeffi¬ 
cients  are  defined  as  [127] 


oq 


d)Q  _  <T)Q 

z  r 


(8) 


where  Q  is  the  molecular  property  in  question,  and  are  the  incident  and 
reflected  fluxes  of  Q,  respectively,  and  is  the  reflected  flux  of  Q  at  complete 
accommodation.  Typical  properties  of  Q  are  tangential  momentum  met ,  normal  mo¬ 
mentum  mcn,  and  kinetic  energy  tangential  mc2/2  and  normal  mc2n/ 2  to  the  surface. 
Here,  m  is  the  molecular  mass,  ct  is  the  tangential  velocity  to  the  surface,  and  cn  is 
the  normal  velocity  to  the  surface.  Accommodation  coefficients  for  these  properties 
are  here  referred  to  as  an,  at,  and  an,  respectively,  so  that 


mct  —  mc[  ct  —  c[ 

(Jt  = 

mct  —  mcw  ct 


(9) 
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where  c't  is  the  post-collisional  tangential  speed,  and  cw  is  the  speed  of  the  wall,  which 
in  this  case  is  stationary.  Similarly, 


(10) 


where  c'n  is  the  post-collisional  normal  speed.  A  useful  relation  that  can  easily  be 
shown  is  at  =  crt(2  —  at). 

This  work  extends  CLL’s  usefulness  by  applying  it  to  adlayer  scattering.  There¬ 
fore,  in  addition  to  the  typical  accommodation  coefficients  needed  by  CLL  ( an  and 
at),  adlayer-specific  coefficients  are  also  required  (an,ad,  &t, ad)- 


2.2.4  Cercignani-Lampis-Lord  Scattering  Kernel. 

The  CLL  scattering  kernel  is  useful  because  it  satisfies  the  principle  of  detailed 
balance  (reciprocity),  it  calculates  sufficiently  realistic  lobular  scattering  angle  and 
energy  distributions,  and  is  well-suited  for  implementation  in  DSMC.  Not  included 
in  this  work  are  Lord’s  extensions  to  CLL,  which  are  diffuse  scattering  with  par¬ 
tial  accommodation  of  translational  kinetic  energy,  vibrational  energy  for  a  rigid- 
rotor/harmonic-oscillator  diatomic  molecule,  partially  diffuse  scattering,  and  vibra¬ 
tional  energy  for  an  anharmonic  oscillator  [98;  99].  Two  of  the  extensions,  which 
consider  vibrational  energy,  are  not  applicable  to  monatomic  xenon.  Future  research 
may  address  these  extensions,  and  apply  them  to  MKS  if  applicable.  Two  system- 
specific  accommodation  coefficients  are  required  for  CLL,  one  for  the  normal  kinetic 
energy  an  and  one  for  the  tangential  kinetic  energy  at.  an  and  at  depend  upon  the 
gas  and  the  surface,  and  are  typically  functions  of  temperature.  Coefficients  have 
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been  measured  experimentally  [52],  as  well  as  predicted  with  MD  simulations  [49].  A 
helpful  explanation  of  the  CLL  kernel  is  given  by  Liou  and  Fang  [96:180-187],  while 
a  slightly  modified  implementation  is  presented  in  [127]. 

The  CLL  scattering  kernel  first  applied  to  DSMC  calculations  by  Lord  is  imple¬ 
mented  in  MKS  [97].  Even  though  CLL  is  applied  here  for  the  first  time  to  adlayer 
scattering,  its  implementation  remains  the  same.  Please  note  that  for  this  section 
only,  all  of  the  velocities  are  normalized  by  the  most  probable  velocity  ^ /2kTs/m , 
consistent  with  Lord’s  development,  where  k  is  the  Boltzmann  constant  and  Ts  is  the 
surface  temperature.  As  a  result,  kinetic  energies  (such  as  u 2)  are  in  units  of  kTs. 
For  each  of  the  tangential  velocity  components  u  and  v,  the  probabilities  that  the 
post-collisional  velocities  fall  between  u'  and  u'  +  du',  and  v'  and  v'  +  dv',  respectively, 
are, 


(11) 


where  P(u  — >  u')  and  P(v  — >  v')  are  the  scattering  kernels.  The  detailed  balance 
relations,  or  reciprocity,  are  then  satisfied  for  both  components, 


exp  {—u2)P(u  — »  v!)  =  exp  {—u'2)P(—u'  — »  —u), 
exp  (— v2)P(y  — >  v')  =  exp  (— v,2)P(— v'  — >  — v ), 


(12) 


The  show  that  the  relations  in  Equation  (12)  hold,  rearrange  the  exponential  terms, 
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and  recall  that  at  =  ot(2  —  at)  =  2 at  —  cp2,  so  that 


exp  (—u2)P(u  — y  u')  = 
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(13) 


and, 
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proving  that 


exp  (— u2)P[u  — >  u')  =  exp  (— u'2)P{— u'  — >  —u) 


(15) 


Similarly, 


exp  ( -v2)P(v  ->  v')  =  exp  (-v'2)P(-v'  -v) 


(16) 


The  normalization  conditions  are  also  satished, 


-OO 

roo 


P(u  — >  u')  du'  =  1, 
P(v  — >■  v')  dv'  =  1, 


(17) 


which  state  that  the  cumulative  probabilities  over  all  scattered  velocities  are  unity. 

For  the  normal  velocity  component  w,  the  probability  that  the  post-collisional 
velocity  is  between  w'  and  w'  +  dw'  is, 


2  ujf 

P(w  — >  w ')  dw'  = -  exp 


w'2  +  (1  -  otn)w2 
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where  P[w  — >  w ')  is  a  scattering  kernel,  and  Iq  is  the  modified  Bessel  function  of  the 
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first  kind.  The  reciprocity  condition  is 


|u?|  exp  (—w2)P(w  — >■  w ')  =  w'  exp  y—wjP  {—w'  — »  — w ) ,  (19) 

where  positive  velocities  are  taken  to  be  away  from  the  surface.  This  condition  is 
different  from  Equation  (12)  since  a  bias  towards  higher  normal  velocities  occurs  in 
the  velocity  distribution  of  impinging  molecules.  The  normalization  condition  is  also 
satisfied, 

P(w  — y  w ')  dw'  =  1.  (20) 

Lord  presented  a  numerical  implementation  of  the  CLL  model  utilizing  a  geomet¬ 
rical  representation,  shown  in  Figure  1  [97].  The  figure  can  be  used  for  both  the 
normal  and  tangential  velocity  components.  Consider  first  the  tangential  velocities  u 
and  v.  Point  P  indicates  the  state  of  the  impinging  molecule  before  it  collides  with 
the  surface.  The  abcissa  and  ordinate  each  represent  one  of  the  tangential  velocities, 
and  so  the  distance  from  the  origin  OR  is  the  tangential  post-collisional  velocity  c!t 
since  R  is  the  state  of  the  reflected  molecule.  P  is  shown  on  the  horizontal  axis 
because  the  velocity  distributions  are  invariant  with  respect  to  rotation.  The  length 
OP  is  |u|.  OQ  is  the  average  value  for  the  post-collisional  speed  (c't)  =  \u\^/l  —  at, 
which  can  be  seen  from  the  fact  that  Equation  (11)  is  a  Gaussian  distribution  with  a 
mean  of  (c't).  The  probability  that  the  final  state  of  the  molecule  lies  in  an  area  dA 
at  (rc,  9C)  is 

P(u  — >  u')P(v  — >  v')  du'  dv'  =  — exp  (  —  —  )  dA ,  (21) 

7TCP  \  at ) 

where  dA  =  rc  drc  d6ci  rc  is  the  distance  QR,  and  9C  is  the  angle  ZPQR.  Note  that 
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the  probability  is  independent  of  9C.  v!  and  v'  are  calculated  from 


6C  =  2tt  Ri  , 


u  =  \u\  \/l  —  at  +  rc  cos  9C, 


v  =  rcsin^c, 


(22) 


where  Ri,  R2  G  [0, 1]  are  uniform  random  numbers.  In  actual  calculations,  one  first 
transforms  velocities  to  the  rotated  reference  frame  in  which  v  =  0.  Once  the  post- 
collisional  velocities  are  obtained,  one  then  transforms  them  back  to  the  original  frame 
of  reference. 

For  the  normal  velocity  component  w,  many  items  stay  the  same.  P  is  the 
molecule’s  initial  state,  the  distance  OR  is  w',  R  is  the  post-collisional  state  of  the 
molecule,  OP  is  |tc|,  OQ  is  |tc| y/1  —  an,  and 


9C  =  2ttR3, 


(23) 


where  R3,  G  [0, 1]  are  uniform  random  numbers,  and  w'  was  determined  with  the 
Law  of  Cosines. 

It  should  be  noted  here  that  MKS  is  not  linked  to  CLL  in  particular,  nor  to 
any  specific  scattering  kernel  nor  method.  MKS  is  modular  in  that  any  scattering 
method  or  kernel  can  be  substituted  in  the  place  of  CLL.  MKS  is  robust  in  that  it 
does  not  need  to  be  modified  to  accommodate  application-specific  scattering  kernels. 
It  has  previously  been  noted  that  CLL  is  not  completely  realistic  [97],  even  though 


Figure  1.  Geometric  representation  of  the  CLL  scattering  model. 

it  has  been  considered  sufficient  for  state-of-the-art  DSMC  calculations  [127].  In 
fact,  other  scattering  kernels  have  been  developed  for  specific  applications  [179;  180]. 
MKS  could  easily  incorporate  other  scattering  kernels  as  long  as  they  are  amenable 
to  DSMC  techniques,  making  it  a  strong 

2.3  Xe-Pt(lll)  System 

A  brief  overview  of  Xe-Pt(lll)  interactions  and  properties  is  included  in  this 
section  since  adsorption  results  will  be  presented  with  respect  to  the  Xe-Pt(lll) 
system.  The  Xe-Pt(lll)  system  includes  attractive  lateral  interactions  among  the 
adatoms  [134],  Xenon  is  a  2D  gas  on  the  surface  if  Ts  is  either  above  the  critical  tem¬ 
perature  Tc,  which  in  this  case  is  120  K,  or  if  Ts  <  Tc  and  0  <  9  <  0.06.  Otherwise, 
the  xenon  phase  is  a  2D  gas  with  either  a  2D  solid  (Ts  <  Ttr)  or  a  2D  liquid  ( Ts  >  Ttr), 
where  Ttr  is  the  triple-point  temperature  [134],  MKS  assumes  no  lateral  adatom  inter¬ 
actions,  and  that  the  adlayer  structure  is  commensurate  up  to  6  =  1. 3  In  reality,  the 
commensurate  monolayer  becomes  saturated  at  6  &  0.4.  The  monolayer  transitions 
through  several  incommensurate  phases  until  the  first  monolayer  is  filled  [74] .  These 

3 A  commensurate  layer  is  one  in  which  adatoms  occupy  energetically-favored  adsorption  sites, 
and  is  dominated  by  the  substrate-adatom  potential.  An  incommensurate  layer,  on  the  other  hand, 
is  determined  by  the  natural  adlayer  interatomic  distance,  or  the  adatom-adatonr  lateral  interac¬ 
tions  [75]. 
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phase  transitions  occur  due  to  competing  forces.  As  long  as  adsorbate-substrate 
interactions  dominate,  the  phase  remains  commensurate.  However,  once  the  lat¬ 
eral  adatom  forces  begin  to  dominate,  the  phase  transitions  to  incommensurate  [74], 
Even  though  one  might  expect  to  see  effects  of  the  phase  transitions  in  experiments, 
none  are  observed  in  Reference  [10].  Arumainayagam  et  al.  caution  that  this  re¬ 
sult  may  be  due  to  the  relatively  high  (1%)  defect  density  on  their  crystal  surface. 
From  Temperature-Programmed  Desorption  (TPD)  spectra,  it  has  been  observed  that 
first-order  desorption  occurs  at  9  <  0.10,  while  for  9  >  0.10,  zero-order  desorption 
is  present  [173].  The  2D  gas  at  low  coverages  explains  the  first-order  behavior,  and 
the  zero-order  desorption  is  explained  by  the  coexisting  2D  solid  and  gas  phases  for 
higher  coverages.  Also,  xenon  atoms  have  a  higher  mobility  than  expected,  due  to 
effective  transfer  of  the  adsorption  energy  into  kinetic  energy  [93]. 

Xe-Pt(lll)  is  a  good  system  with  which  to  initially  compare  an  adsorption  method 
because  it  is  relatively  simple.  Xenon  is  an  inert  gas,  and  the  Pt(lll)  samples  may  be 
single-crystal,  clean,  and  well-defined  metal  surfaces  with  a  relatively  smooth  PES. 
The  adsorption  of  xenon  on  Pt(lll)  is  a  non-activated  process,  meaning  that  no 
chemical  reaction  occurs  during  adsorption.  It  is  a  weakly-interacting  gas-surface 
system  since  the  energy  is  transferred  to  the  surface  by  exciting  phonons.  No  internal 
energy  transfer  occurs  in  the  incident  species. 

2.3.1  Xe-Pt(lll)  Potentials. 

In  simulating  the  dynamics  of  the  Xe-Pt(lll)  system,  different  potentials  have 
been  investigated,  namely  the  LJ,  the  Barker-Rettner  (BR),  and  the  Morse  potentials. 
The  performance  of  MKS  will  be  compared  with  calculations  from  the  BR  and  Morse 
potentials  in  Section  2.4.2,  and  since  the  Morse  potential  is  based  on  the  LJ  potential, 
all  three  potentials  will  be  briefly  discussed  here. 
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The  LJ  potential  Vlj  has  two  adjustable  parameters  (eLJ  and  <rLJ),  and  is  of  the 


form 


VLj  =  4£l 


J  LJ 

r 


J  LJ 

r 


(24) 


where  <tlj  is  the  intermolecular  distance  at  which  the  potential  is  zero,  and  r  is  the 
intermolecular  distance. 

The  Morse  potential  Vm  has  three  adjustable  parameters  (eM,  crM,  and  r0),  and  is 
of  the  form  [9] 


VM  =  £m  {exp  [— 2crM(r  -  r0)]  -  2  exp  [~aM(r  -  r0)]}  ,  (25) 


where  eM  is  the  well-depth  parameter,  crM  is  a  parameter  which  adjusts  the  steepness 
of  the  potential,  and  vq  is  the  equilibrium  intermolecular  distance. 

The  BR  potential  has  a  more  complicated  form  than  the  LJ  and  Morse  potentials. 
It  was  developed  specifically  for  the  Xe-Pt(lll)  system,  and  is  intended  for  MD 
simulations.  The  BR  potential  RBr  has  nine  adjustable  parameters  aq,  c6,  hs, 
rq,  Ag,  W,  a,  and  5)  with  the  following  form  [13;  172]: 

VBR  =  £  Mr,)  +  „(r')l  +  (26) 


where 


u{ri)  =  A\  exp(— air*), 


(27) 


and 


v(r[)  = 


c6  G(r') 
„/6 


(28) 


where  Ai,  oq,  and  cq  are  adjustable  parameters,  rq  is  the  distance  between  the  gas 
molecule  and  the  ith  bulk  molecule,  and  r[  is  centered  a  distance  hs  above  the  bulk 
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molecule  with 


(29) 


r't  =  (x,  -  Xi)1  +  (ya  -  ft  y  +  (zg  -Zi-  h,)‘, 


where  (ay,  yi:  )  is  the  position  of  the  bulk  molecule,  the  coordinates  of  the  gas 
molecule  are  ( xg ,  yg,  zg),  and  hs  is  an  adjustable  height  parameter.  G(r )  is  a  damping 
function  defined  as 


G(r)  = 


exp 


^-1 

r 


for  r  >  ry, 
for  r  <  ry, 


(30) 


where  ?y  is  an  adjustable  parameter.  The  term  V (zfre)  is  defined  as 


Aq  W  exp  (-c<ve) 
V(z™e)  =  9  V  9  ’ 


W  +  Ag  exp  (— az^ve)  ’ 


(31) 


where  Ag,  W,  and  a  are  adjustable  parameters,  and  z*ve  is  the  height  above  the 


‘local- average”  surface 


ave  _  _  ave 

^ g  —  ^ 9  * gs  5 


(32) 


where  z™e  is  a  weighted  average  over  the  surface  molecules, 


„ave  =  Es  <l>(ri) 
Zgs  — 


Zttr.)  ’ 


(33) 


with  the  subscript  s  indicating  that  the  summation  in  the  numerator  is  performed 
only  for  the  surface  molecules.  <j)(r)  is  a  weighting  function  given  by 


(j){r)  =  exp  (— <5r2)  , 


(34) 


where  5  is  an  adjustable  parameter. 
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2.4  Results 


To  verify  that  MKS  may  be  utilized  to  simulate  adsorption  and  scattering  in 
DSMC  simulations,  it  will  be  shown  in  this  section  that  (1)  not  only  is  the  integrity 
of  CLL  maintained,  but  its  performance  is  actually  improved,  (2)  initial  adsorp¬ 
tion  probabilities  S0  compare  favorably  with  experiment,  and  (3)  coverage-dependent 
probabilities  agree  well  with  data. 

2.4.1  Scattering. 

The  scattering  properties  of  CLL  are  improved  when  incorporated  into  MKS. 
Figures  2  and  3  compare  experimental  data  with  MKS  and  CLL  calculations  for 
scattering  angles  and  E'tr  [13] .  The  incident  angle  is  45°  and  the  incident  translational 
energy  is  48.2  kJ  mol-1.  The  CLL  and  MKS  results  were  sampled  from  10'  simulated 
particles.  Only  molecules  scattered  in-plane  (±1°)  are  reported,  to  match  what  was 
measured  by  experiment  [143;  147;  148].  In  Figure  2,  the  in-plane  scattering  angle 
is  measured  from  the  surface  normal.  The  scattered  intensity  for  MKS  has  a  tighter 
distribution  than  CLL,  and  approaches  experimental  values.  In  Figure  3,  E[r  is  the 
average  total  energy  in  that  scattering  direction.  At  higher  angles,  MKS  approaches 
CLL,  while  at  lower  angles,  MKS  follows  the  upward  trend  seen  from  experiment.  It 
is  interesting  to  note  that  when  the  gas  molecule  collides  with  the  surface,  it  may 
scatter  such  that  E'tr  >  EtI.  This  phenomenon  is  termed  super-elastic  scattering,  and 
describes  the  case  where  the  gas  molecule  picks  up  thermal  energy  from  the  surface 
during  the  collision. 

It  would  be  preferable  to  compare  scattering  data  at  the  same  surface  temperature 
used  in  Sections  2.4.2  and  2.4.3  (i.e. ,  95  K).  However,  there  is  little  scattering  data 
for  the  Xe-Pt(lll)  system  in  the  literature.  In  Reference  [13],  scattering  data  were 
reported  for  Xe-Pt(lll)  at  Ts  =  800  K,  which  is  much  higher  than  Ts  =  95  K.  In 
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Figure  2.  Scattered  intensity  versus  in-plane  scattering  angle  for  Xe-Pt(lll).  Calcula¬ 
tions  from  CLL  and  MKS  are  compared  with  experiment.  The  surface  temperature  is 
Ts  —  800  K.  The  angle  of  incidence  6i  and  the  in-plane  scattering  angle  are  measured 
from  the  surface  normal.  The  angle  of  incidence  Qi  is  45°  with  an  incident  energy 
Etr  —  48.2  kJ  mol-1.  The  vertical  line  indicates  purely  specular  scattering.  CLL  and 
MKS  simulations  were  performed  with  an  =  0.74,  att  —  0.60,  eLJ  =  2.142  kJ  mol-1,  and 
-Eads  =  25.9  kJ  mol-1.  — ,  MKS; - ,  CLL;  o,  experimental  data  from  [13],  Figure  7. 
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Figure  3.  Final  translational  energy  E'tr  versus  in-plane  scattering  angle  for  Xe-Pt(lll). 
Calculations  from  CLL  and  MKS  are  compared  with  experiment.  The  surface  temper¬ 
ature  is  Ts  —  800  K.  The  angle  of  incidence  0i  is  45°.  The  horizontal  line  represents 
the  incident  energy  Etr  —  48.2  kJ  mol-1.  CLL  and  MKS  simulations  were  performed 
with  OLn  —  0.74,  at  —  0.60,  eLJ  =  2.142  kJ  mol-1,  and  Ea ds  =  25.9  kJ  mol-1.  — ,  MKS; 
- ,  CLL;  o,  experimental  data  from  [13],  Figure  7. 


spite  of  this  deficiency,  comparison  with  experiment  at  the  elevated  temperature  still 
provides  an  opportunity  to  consider  how  MKS  performs. 


2.4.2  Initial  Adsorption  Probabilities. 

Arumainayagam  et  al.  reported  experimental  initial  adsorption  probabilities  Sq 
for  the  Xe-Pt(lll)  system  using  supersonic  atomic  beam  techniques  [9].  The  surface 
temperature  of  Pt(lll)  was  held  at  Ts  =  95  K  while  the  incident  translational  energy 
Etr  of  the  xenon  atoms  varied  from  7  to  63  kJ  mol-1  at  incident  angles  from  0° 
to  60°  with  respect  to  the  surface  normal.  At  these  incident  energies,  the  radius 
parameter  1Z  is  between  1.6  and  1.7  (2.52  A  <  RcS  <  2.65  A,  d0  =  2.70  A,  and 
Rc  =  1.56  A),  indicating  that  this  system  falls  under  the  thermal  scattering  regime 
and  that  MKS  is  applicable.  Sq  was  measured  within  an  experimental  error  of  ±0.02, 
and  the  coverage  6  noise  level  was  about  0.01  ML. 
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MKS  simulations  are  compared  with  data  in  Figure  4  [9].  Each  point  on  the 
MKS  curve  is  from  a  sample  of  105  molecules  impinging  on  the  surface  with  incident 
normal  and  in-plane  tangential  velocities  cn  =  ccosOi  and  ct  =  csin  0j,  respectively, 
where  c  =  \J2 Etr/m  is  the  total  incident  velocity,  and  m  is  the  gas  molecule  molar 
mass.  Points  were  chosen  at  11  different  angles  equally  spaced  from  0°  to  90°,  and 
11  different  specific  energies  equally  spaced  from  0  to  63  kJ  mol-1.  For  nm  =  2.0,  all 
of  the  points  fell  onto  a  single  curve.  The  values  for  an  and  at  were  chosen  to  best 
match  experiment  in  Figure  6,  and  are  both  within  the  published  ranges  of  0.7  —  1.0 
and  0.6  —  0.8,  respectively  [49].  The  simulations  were  less  sensitive  to  at  than  to  an, 
indicating  that  the  MKS  method  would  perhaps  need  to  be  modified  to  handle  larger 
surface  corrugation;  the  PES  for  Xe-Pt(lll)  is  relatively  smooth.  anjad  and  aya d  were 
not  specified,  since  the  surface  was  assumed  to  be  clean.  The  other  parameters  were 
taken  from  the  literature  [9]:  eLJ  =  2.142  kJ  mol-1,  and  Ea ds  =  25.9  kJ  mol-1,  where 
f7ads  is  the  molar  energy  of  adsorption.  The  surface  temperature  was  Ts  =  95  K.  The 
exponent  nm  from  Etr  cosnm  6i  is  not  based  in  theory,  but  has  been  found  to  place 
So  data  points  on  a  single  curve.  nm  =  2.0  was  chosen  in  Figure  4  to  highlight  the 
MKS  curve.  In  Figure  5,  MKS  calculations  are  compared  with  values  obtained  by 
Weaver  et  al.  [172]  from  Morse  and  BR  potentials,  with  nm  =  1.8,  105  simulated 
impinging  molecules  per  point  at  50  different  angles  equally  spaced  from  0°  to  90°, 
and  50  different  specific  energies  equally  spaced  from  0  to  63  kJ  mol-1,  with  the  other 
parameters  the  same  as  in  Figure  4.  The  BR  and  Morse  potential  parameters  are 
listed  in  the  caption  to  Figure  5.  Admittedly,  MKS  does  not  recreate  S0  well,  even 
though  it  will  be  shown  in  Figure  6  to  match  S(9)  much  better.  For  the  results  in 
Figures  4  and  5,  adsorption  is  occurring  as  a  2D  gas. 
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Figure  4.  Initial  adsorption  probability  So  versus  modified  energy  Etrcos20i  for  Xe- 
Pt(lll).  The  surface  temperature  is  Ts  =  95  K.  The  curve  is  from  MKS  simulations 
with  an  —  0.74,  at  —  0.60,  eLJ  =  2.142  kJ  mol-1,  and  Ea ds  =  25.9  kJ  mol-1.  — ,  MKS; 
A,  7  kJ  mol-1  (exp);  □,  18  kJ  mol-1  (exp);  o,  37  kJ  mol-1  (exp);  V,  63  kJ  mol-1  (exp). 
Experimental  data  are  from  [9],  Figure  4. 


2.4.3  Coverage-Dependent  Adsorption  Probabilities. 

The  adsorption  probability  S(9)  is  also  a  function  of  coverage.  Adsorption  proba¬ 
bility  versus  coverage  for  five  different  initial  conditions  is  shown  in  Figure  6;  the  data 
are  from  [10]  with  uncertainties  of  ±3  K  for  Ts,  and  ±0.03  for  S(9).  The  parameters 
are  the  same  as  those  from  Figure  4.  In  addition,  an^  =  0.89  and  cpiad  =  0.60  were 
chosen  in  order  to  best  fit  MKS  results  to  experiment.  S{9)  was  not  sensitive  to  the 
choice  of  a^ad-  S{9)  increases  with  coverage,  requiring  anjad  >  «n.  This  finding  is 
consistent  with  observations  made  in  [10]  where  the  adsorption  probability  on  covered 
surfaces  was  greater  than  the  adsorption  probability  on  a  clean  surface.  The  results 
are  a  2D  gas  where  9  <  0.02,  and  a  2D  solid  combined  with  a  2D  gas  otherwise. 

A  sensitivity  analysis  was  performed  to  determine  the  optimal  value  for  the  post- 
collisional  energy  threshold  value,  and  it  was  verified  in  this  work  that  2eLJ  provides 
satisfactory  results.  Yamanishi  et  al.  claim  that  when  a  gas  molecule  scatters  off  of 
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Figure  5.  Initial  adsorption  probability  So  versus  modified  energy  Etrcos186i  for  Xe- 
Pt(lll).  The  surface  temperature  is  Ts  =  95  K.  The  curve  is  from  MKS  simulations 
with  an  —  0.74,  a.t  =  0.60,  eLJ  =  2.142  kJ  mol-1,  and  Ea ds  =  25.9  kJ  mol-1.  *,  MKS; 
•  ,  BR;  O,  Morse.  Morse  potential  parameters  are  [172]:  £M  =  2.628  kJ  mol-1,  <rM  = 
1.05  A  1,  and  ro  =  3.20  A.  BR  parameters  are  [13]:  Ag  —  3.084  X  106  kJ  mol-1,  a  — 
4.25  A-1,  W  —  96  kJ  mol-1,  Ax  =  5.718  X  1011  kJ  mol-1,  a±  =  9.00  A-1,  S  =  0.22  A-2, 
Cq  =  7.981  X  103  kJ  A6  mol- 1,  hs  —  1.83  A,  and  ri  =  6.8247  A  [14].  BR  and  Morse 
calculations  are  from  [172],  Figure  3. 
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6i  =  0  °  6i  =  0  °  6i  =  30  ° 


0i  =  0  °  6i  =  45  ° 


Figure  6.  Adsorption  probability  S(6)  versus  coverage  6  for  Xe-Pt(lll)  at  different 
incident  kinetic  energies  and  angles.  The  surface  temperature  is  Ts  =  95  K.  6  is 
measured  in  monolayers.  MKS  simulations  were  performed  with  an  —  0.74,  at  =  0.60, 
c^n, ad  =  0.89,  at:& d  =  0.60,  £LJ  —  2.142  kJ  mol-1,  and  Ea ds  =  25.9  kJ  mol-1.  — ,  MKS; 
o,  experimental  data  is  from  [10],  Figure  1. 


39 


a  surface,  it  has  an  average  potential  energy  of  2eLJ,  and  so  2eLJ  was  the  threshold 
value  that  they  chose  [180].  For  the  MKS  method,  values  higher  than  2eLJ  were 
discovered  to  improve  post-collisional  energies  and  scattering  angles  when  compared 
with  experiment,  but  caused  the  adsorption  probabilities  to  deviate  even  more  from 
experimental  values.  For  a  value  of  5eLJ,  calculated  post-collisional  energies  and 
scattering  angles  agreed  well  with  experiment,  while  adsorption  probabilities  increased 
dramatically  and  were  too  large  compared  with  experimental  values.  On  the  other 
hand,  threshold  values  less  than  2eLJ  did  not  improve  any  of  the  MKS  results. 

2.5  Physisorption  Summary 

A  new  method,  the  Modified  Kisliuk  with  Scattering  method,  has  successfully 
been  developed  for  calculating  adsorption  probabilities  and  scattering  properties  in 
Direct  Simulation  Monte  Carlo  applications.  The  Modified  Kisliuk  with  Scattering 
method  incorporates  a  commonly-employed  scattering  kernel,  the  Cercignani-Lampis- 
Lord  kernel,  with  the  underlying  assumptions  of  an  existing  adsorption  model,  the 
Modified  Kisliuk  model.  The  results  show  that,  for  the  Xe-Pt(lll)  system,  the  new 
method  was  not  only  able  to  reproduce  adsorption  probabilities  as  a  function  of 
coverage,  but  it  also  improved  the  scattering  properties  of  the  Cercignani-Lampis- 
Lord  kernel. 

This  work  provides  two  significant  contributions  to  the  Cercignani-Lampis-Lord 
kernel.  First,  the  kernel  was  applied  to  scattering  off  of  an  adlayer  with  favorable 
results.  Previously,  work  with  this  kernel  has  been  confined  to  the  gas  scattering 
off  of  the  primary  surface.  This  is  the  first  time  that  gas-adlayer  interactions  have 
been  modeled  with  the  Cercignani-Lampis-Lord  kernel.  Second,  optimal  values  for 
the  accommodation  coefficients  were  found  by  comparing  results  with  experiment. 
The  Modified  Kisliuk  with  Scattering  method  thus  provides  an  additional  method  for 
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determining  appropriate  accommodation  coefficients. 

Direct  Simulation  Monte  Carlo  techniques  have  now  been  shown  to  be  capable 
of  simulating  surface  effects,  accomplished  by  successfully  modeling  adsorption  and 
scattering  for  a  physisorptive  system.  The  next  step,  examined  in  Chapter  III,  will 
be  to  demonstrate  that  desorption  in  a  chemisorptive  system  can  also  be  modeled  for 
Direct  Simulation  Monte  Carlo. 
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III.  Chemisorption 


Chemisorption  is  much  more  complex  than  physisorption,  and  will  therefore  re¬ 
quire  different  modeling  methods  than  those  presented  in  Chapter  II.  In  addition, 
Chapter  II  was  concerned  with  adsorption,  whereas  desorption  will  now  be  exam¬ 
ined,  the  microscopically  inverse  process  of  adsorption.  However,  desorption  is  more 
complicated  than  adsorption  for  at  least  two  reasons.  First,  the  timing  of  desorption 
events  must  be  determined  accurately  for  a  temporal  simulation.  Adsorption  events 
simply  occur  when  a  gas  molecule  impinges  the  surface,  and  in  a  DSMC  code,  the 
positions  and  velocities  of  all  the  gas  molecules  are  known  and  tracked.  Therefore,  the 
timing  of  adsorption  events  is  not  an  issue.  However,  determining  when  an  admolecule 
will  (statistically)  attain  a  sufficiently  high  enough  energy  to  desorb,  obtained  from 
the  surface  thermal  energy,  is  not  trivial.  Second,  there  are  many  desorption  path¬ 
ways  for  the  same  adsorbate.  As  discussed  in  Chapter  I,  desorption  may  be  caused 
by  surface  thermal  energy  (LH),  strong  EM  fields  (FID),  atom/molecule  (ER),  ion 
(IID),  photon  (PID),  or  electron  (ESD)  surface  impact,  and  each  of  these  pathways 
involves  separate  physical  phenomena. 

Desorption  induced  by  surface  thermal  energy  (LH)  is  thought  to  be  a  signifi¬ 
cant  contributor  to  hydrogen-gas  contamination  in  magnetron  devices,  therefore  that 
pathway  is  investigated  here  for  the  relevant  H2-Cu(100)  system  (Section  3.1)  in 
the  following  manner.  Each  desorbing  molecule  is  modeled  with  the  CT  method 
(Section  3.3)  as  it  is  influenced  by  a  6D  PES  (Section  3.2).  Initial  conditions  are 
assigned  at  the  average  TS  locations  (Section  3.6)  according  to  a  truncated  Maxwell- 
Boltzmann  (MB)  distribution  (Section  3.7).  Then,  the  molecule  is  simulated  forwards 
and  backwards  in  time  according  to  Keck’s  method  (Section  3.5)  to  determine  if  that 
specific  trajectory  represents  a  valid  desorption  event.  Since  simulations  are  per¬ 
formed  at  the  molecular  scale,  all  of  the  relevant  equations  are  non-dimensionalized 
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(Section  3.4).  The  timing  of  desorption  events  is  determined  via  a  Poisson  process 
given  a  desorption  rate  (Section  3.8).  Outputs  of  this  method  include  the  desorption 
angle,  the  translational  energy,  and  the  internal  energies  of  the  desorbing  H2  molecule 
(Section  3.9). 

This  work  is  a  significant  step  forward  for  at  least  seven  reasons.  First,  desorption 
modeling  has  not  been  available  in  DSMC  applications  until  now.  Desorption  has  been 
ignored  in  the  DSMC  community  due  to  its  complexity.  The  methods  presented  here 
can  be  directly  applied  to  DSMC  applications  as  a  desorption  boundary  condition. 
Second,  the  timing  of  desorption  events  can  now  be  coupled  with  CT  simulations 
to  model  temporal  desorption,  a  necessary  requirement  for  system-level  desorption 
simulations.  Third,  a  non-dimensionalization  scheme  is  developed  for  use  with  the 
CT  method.  Not  a  single  non-dimensionalization  scheme  was  found  in  the  literature. 
Fourth,  this  work  considers  all  of  the  contributing  TS  locations  and  weights  them  ac¬ 
cordingly.  Previous  work  only  considered  one  TS  location  when  implementing  Keck’s 
method.  Fifth,  a  truncated  MB  distribution  is  developed  (with  its  accompanying 
accept-reject  form)  to  determine  initial  conditions  at  the  TS  locations.  Sixth,  state- 
resolved  average  translational  energies  are  accurately  modeled  without  the  need  for 
resource-intensive  Gaussian  Weighting  (GW).  Finally,  the  absorption  energy  barrier 
is  shown  to  significantly  contribute  to  the  translational  energy  of  desorbing  molecules, 
but  has  little  effect  on  their  rotational  and  vibrational  energies. 

3.1  H2-Cu(100)  System 

H2-CU  is  the  simplest  chemisorption  desorption  system  to  model  [126],  and  the 
most  thoroughly  studied  example  of  activated  adsorption  [114],  It  is  a  fully-activated 
system,  with  activation  barriers  ranging  from  0.5  to  1.0  eV  [11;  126].  The  strength 
of  these  barriers  ensures  that  there  is  virtually  no  desorption  (adsorption)  on  H2-CU 
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Figure  7.  The  unit  cube  for  an  FCC  crystal  structure.  The  lattice  parameter  a  is 
defined  as  the  average  equilibrium  length  of  the  unit  cube.  The  top  face  is  highlighted 
to  show  the  (100)  surface. 

at  room  temperature  [5].  Since  these  average  barriers  have  been  shown  to  be  only 
weakly  dependent  on  surface  coverage  [166],  the  assumption  that  desorption  events 
occur  independently  will  be  asserted  here.  This  assumption  equivalently  describes  a 
system  in  which  desorption  events  are  dependent  upon  one  another,  but  where  there 
is  (almost)  zero  surface  coverage. 

When  describing  the  surface  of  a  crystal,  it  is  important  to  identify  how  the 
crystal  is  cut  in  reference  to  the  unit  cube.  Copper  exhibits  a  face-centered  cubic 
(FCC)  crystalline  structure,  shown  in  Figure  7  as  a  unit  cube.  Copper  atoms  are 
placed  at  each  corner  of  the  cube  and  at  the  center  of  each  of  the  six  faces.  As  a  brief 
introduction  to  the  nomenclature,  the  (111)  surface  is  found  by  cutting  the  unit  cube 
parallel  to  the  plane  defined  by  the  vector  ( x,y,z )  =  (1,1,1)  and  passing  through 
the  points  (x,  y,  z )  =  (a,  0,  0),  (0,  a,  0),  and  (0,  0,  a),  where  a  is  the  lattice  parameter, 
or  the  average  equilibrium  length  of  the  unit  cube.  Similarly,  the  (100)  face  (see 
Figure  7)  is  parallel  to  the  plane  defined  by  the  vector  (1,  0,  0)  and  the  point  (1,  0,  0). 
Note  that  since  all  six  faces  of  an  FCC  crystal  are  identical,  the  (100),  (010),  and 
(001)  surfaces  are  equivalent. 
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Figure  8.  A  comparison  of  the  unit  cube  and  cuboid  on  the  Cu(100)  surface.  Filled 
circles  represent  the  copper  atoms  in  the  top  layer,  while  empty  circles  represent  the 
layer  of  copper  atoms  immediately  below  the  surface  layer.  The  unit  cube  from  Figure  7 
is  shown  with  dashed  lines,  including  the  typcial  lattice  parameter  a.  However,  the 
smaller  unit  cuboid  enclosed  within  the  solid  lines,  with  its  corresponding  distance  <n, 
is  the  unit  cell  to  which  reference  is  made  in  this  work. 


The  unit  surface  cell  utilized  in  this  work,  shown  in  Figure  8,  differs  slightly  from 
the  unit  cube  defined  in  Figure  7.  The  unit  cell  is  a  rectangular  cuboid,  the  smallest 
repeating  unit  on  the  (100)  FCC  surface.  By  inspection,  the  length  ep  is  calculated 
from  the  relation  cq  =  a/ y/2.  The  experimental  value  a  for  copper  is  a  =  3.609  A  [174], 
so  that  ai  =  2.552  A.1 


3.2  Potential  Energy  Surface  (PES) 

The  interaction  between  a  molecule  and  a  surface  is  described  by  the  PES,  from 
which  one  can  calculate  the  relevant  forces.  Then,  with  appropriate  initial  conditions, 
the  molecule’s  position  and  velocity  can  also  be  determined.  Diatomic  hydrogen  has 
six  degrees  of  freedom,  and  therefore  a  6D  PES  is  required  for  a  full  description  of 
the  molecular  motion  [44;  45]. 

Wiesenekker  et  al.  [174]  developed  a  6D  PES  to  describe  the  molecule-surface 

1In  units  of  the  Bohr  radius  a o,  a  =  6.820  a o  and  a;  =  4.822  oq. 
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Figure  9.  An  example  2D  slice  of  the  full  PES.  Taken  at  x  —  at,  y  =  ai/ 2,  6  =  90°, 
and  <f)  —  90°,  copresponding  to  the  case  where  hydrogen  atoms  are  adsorbed  on  top  of 
adjacent  copper  atoms,  and  desorbed  as  H2  above  a  bridge  site.  The  TS  is  located  at  the 
saddle  point,  indicated  by  an  X.  The  contours  represent  potential  energy  in  units  of  eV. 
The  reaction  path  would  follow  the  path  of  lowest  potential  energy  from  the  adsorbed 
state  (r  =  ai )  to  the  desorbed  state  ( z  =  1.87  ai ),  passing  through  the  TS.  Note  that 
the  H2-Cu(100)  system  exhibits  early  (late)  barriers  to  desorption  (adsorption)  [139]. 


interaction  for  the  H2-Cu(100)  system.  The  PES  is  based  on  slab  calculations  with 
the  GGA  in  DFT  with  an  RS  assumption  (bulk  atoms  are  frozen  in  time  and  space), 
and  is  expressed  as  an  analytical  fit,  allowing  for  a  high  degree  of  flexibility,  and 
bypassing  the  need  for  computationally  intensive  on-the-fly  DFT  calculations.  The 
full  6D  PES  is  achieved  by  splicing  together  eight  2D  PES  slices  [175;  176],  which  are 
individually  anchored  to  experimental  results.  This  6D  PES  was  specifically  tailored 
to  be  accurate  in  the  region  between  the  gas  phase  and  the  TS  (the  location  or  zone 
where  the  individual  adatoms  combine  into  a  molecule).  Within  the  TS,  the  PES  is 
fairly  accurate,  and  between  the  TS  and  the  surface,  the  PES  is  inaccurate.  However, 
since  results  will  focus  on  desorption  properties  (with  initial  conditions  specified  at 
the  TS),  the  region  where  the  PES  is  inaccurate  will  not  play  a  large  enough  role  to 
be  a  concern. 

The  6D  Wiesenekker  PES  is  constructed  in  two  stages.  First,  eight  2D  PES  slices 
are  calculated  as  functions  of  r  (interatomic  bond  length)  and  z  (shortest  distance 
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from  the  center  of  mass  of  the  molecule  to  the  surface)  for  various  molecular  orien¬ 
tations  and  desorption  paths.  The  eight  2D  PESs  are  then  combined  into  a  6D  PES 
by  solving  sets  of  linear  equations.  The  details  of  this  process  are  provided  below  for 
the  convenience  of  the  reader. 

Four  of  the  eight  2D  PES  slices  (V^go,  14tg0,  VAgo,  and  V)b9o)  are  created  by 
summing  the  contributions  from  a  two-body  potential  V2  and  a  three-body  potential 


V (r,  z)  =  V2(r,  z)  +  V3(r,  z), 


(35) 


where  the  two-  and  three-body  potentials  are  constructed  similarly, 


C<Co-AC, 


v*(r>  z )  /c(C )V2A(r,  z)  +  [1  -  fe(Q]V2B(r,  z),  Co  -  AC  <  C  <  Co  +  AC,  (36) 


C  >  Co  +  AC, 


and 


U(r.  z)  =  ^  /c(C)iy (r,  z)  +  (1  -  fc(Q]V3B(r,  z),  Co  -  A(  <  C  <  Co  +  AC,  (37) 


and  fc( C)  is  the  switching  function,  defined  by 


(38) 


with 


7T  rc-(Co-AC) 
2  [  AC 


(39) 
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and 

C— (i^).  (40) 

Values  for  the  above  parameters  are  given  in  Reference  [174]  as  zre f  =  18.3  a0,  rre f  = 
11.0  a0,  Co  =  61.5°,  and  AC  =  2.5°.  The  two-body  expressions  VA  and  V2B  are 
superpositions  of  attractive  (Ktt)  and  repulsive  (KeP)  potentials, 


V2A(r,z ) 

=  Ktt  (r)  +  Kep  (z), 

(41) 

V2 B(r,z) 

—  2  Ktt  (z)  +  Kep  (r)- 

(42) 

V2a  represents  the  H2  molecule  as  it  desorbs,  so  that  the  molecule  experiences  an 
attractive  force  between  the  two  hydrogen  atoms  but  a  repulsive  interaction  with  the 
surface.  On  the  other  hand,  V2  represents  the  two  H  atoms  before  they  combine 
into  an  H2  molecule,  so  that  both  an  attractive  interaction  with  the  surface,  and  a 
repulsive  interaction  between  the  atoms,  exist.  The  attractive  potential  is  fit  to  a 
modified  Rydberg  form, 

Ktt (p)  =  -De  (1.0  +  diPatt  +  a2pltt  +  a3Patt)  exp(— a4patt),  (43) 

where  patt  =  r  —  re  in  Equation  (41),  but  patt  =  z  —  ze  in  Equation  (42),  and  the  Pauli 
repulsion  Vrep  is  expressed  as 


Kep(Prep)  ®  exp(  bpIe p),  (4 4) 

where  prep  =  r  in  Equation  (41),  prep  =  z  in  Equation  (42),  and  the  remaining 
parameters  are  found  in  Tables  II  and  III  in  Reference  [174].  The  actual  fitting  is 
contained  within  the  three-body  expressions  VA  and  VB,  which  represent  the  energy 
differences  between  the  GGA  energy  and  the  two-body  potential,  fitted  with  a  least- 
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squares  procedure  in  the  following  expressions, 


V3A(r,  z)  =  P(su  s2)  [1.0  -  tanh(7isi)]  [1.0  -  tanh(72s2)] ,  (45) 

Vf(r,  z)  =  P(su  s2)  [1.0  -  tanh(7isi)]  [1.0  -  tanh(72s2)] ,  (46) 

where  P(s3,s2)  is  the  analytical  curve  fit, 

P(S  1,  S2)  =C0  +  CiSi  +  c2s2 

+  Ci\s\  +  Ci2SiS2  +  C22S^ 

+  Cmsf  +  Cn2SfS2  +  Ci22SiS2  +  C222S2  (^) 

+  CnnSi  +  Cm2sfs2  +  Cn22SfS2  +  Ci222SiS2  +  C2222s(> 

+  Cnii2S^S2  +  Cm22sfs2  +  Cn222SfS2  +  Ci2222SiS2, 

with  si  =  r  —  r0  and  s2  =  z  —  Zq.  Note  that  if  r  <  r3b,  then  r  =  r3b.  Likewise,  if  <  £3b, 

then  z  =  z3b-  All  of  the  parameters  for  V  ^  and  V  ^  are  taken  from  Tables  IV  and  V 

in  Reference  [174],  respectively. 

The  remaining  four  2D  PES  slices  (Khuo,  Vbtuo,  V/(,6i40,  and  V*bi4o )  are  also  super¬ 
positions  of  two-body  (V2)  and  three-body  (V3)  interactions, 

V  =  V2(r,z,0)  +  V3(r,z).  (48) 


The  two-body  potential  V2  is 


V2(r,  z,  6)  =  I4tt (r )  +  Viep(z,  9),  (49) 

where  Vatt  is  taken  from  Equation  (43)  with  p  =  r  —  re  and  the  parameters  from  the 
column  labeled  H2  in  Table  II,  Reference  [174].  The  repulsive  potential  is  dependent 
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upon  9, 


KeP (z,  0)  =  ^  [exp(-fezi)  +  exp(— 5z2)]  ,  (50) 

where  z3  and  z2  are  the  ^-coordinates  (height  above  the  surface)  of  the  individual 
atoms, 


Zi  =  z  — 


r  cos  6 
2 


z2  —  z  + 


r  cos  6 
2 


(51) 


and  the  parameters  a  and  b  are  found  in  the  column  for  Cu-H2  interactions  in  Ta¬ 
ble  Ill,  Reference  [174],  The  three-body  potential  VA  is  again  calculated  from  Equa¬ 
tions  (45)  and  (47),  with  Si  =  r  —  r0  and  s2  =  z  —  z0.  If  r  <  r3b  then  r  =  r3b,  and 
if  z  <  z3b  then  z  =  z3b  (parameters  taken  from  Table  VI  in  Reference  [174]).  V3  is 
determined  from, 


V3(r,z)  =  < 


V3A(r,z ), 
m)V3A(r,z), 
0, 


f  <  £o  -  Af , 

&  -  <  e  <  e0  +  Af , 

e  >  &  +  Ae, 


(52) 


with  the  help  of  another  switching  function  fd(£),  similar  to  fc( C),  where 


m) 


1  +  COS  X2 
2 


(53) 


7T  re-fco-AQ 

2  [  A£ 


(54) 


and 

£  =  tan-1  ( - — .  (55) 

\r  —  r  ref)2  / 

Values  for  the  above  parameters  are  given  as  zre fi2  =  14.0  ao,  rrefi2  =  15.0  ao,  £o  = 
46.15°,  and  A£  =  1.15°  in  Reference  [174].  As  a  side  note,  there  is  an  error  in  Figure  5a 
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in  Wiesenekker  et  al.  [174].  That  plot  does  not  represent  calculations  with  the  pa¬ 
rameter  z3b,  as  it  should.  Figures  3a,  3b,  4a,  4b,  5a,  and  5b  in  Wiesenekker  et  al.  [174] 
can  be  recreated  by  plotting  Vbh90  (r,  z),  Vbt90  (r,  z),  Vhb90  (r,  z),  Vmo  (r,  z),  VbhU0  (r,  2,  9), 
and  Vbtuo(r,  z,  6),  respectively. 

The  2D  expansion  coefficients  are  constructed  from  the  various  2D  PES  slices. 
They  are  grouped  according  to  whether  desorption  occurs  over  a  bridge  site, 


V20 b(r,  z,  9)  = 


Vqo b(r,  z)  = 


Vfefci4o(y,  Z,  9)  +  Vhtuojr,  z,  6)  -  Vbh90 (r,  z)  -  Vbm (r,  z) 
2  [y°(0  =  140.8°)  -  Y$(6  =  90°)] 

Vbh90 (r,  z)  +  14*90(7-,  z)  -  V20 b(r,  z)Y£(9  =  90°) 

2^0° 


V2eb(r,  z)  = 


Vbhdojr,  z)  -  14*90(7-,  z) 

2 Y2e(6  =  90°,  0  =  0°)  ’ 


(56) 


a  hollow  site, 


or  a  top  site, 


V20h{r,z,9)  = 


V00h(r,z)  = 


VhbU0(r,  z,  9)  -  Vh 690  (r,  z) 


V20  t(r,z,9)  = 


Vqo  t(r,z)  = 


Y2y  = 

140.8°) 

O 

O 

1 

Vhb9o(r , 

2) 

-V2Oh{r,z)Y2°(9  =  90°) 

v° 

14,140(7-, 

z,9) 

-  14690(7-,  z) 

Y?(0  = 

:  140.8°) 

1 

QO 

O 

O 

Vmo  (t1  , 

*)■ 

-V20t(r,  z)Y${9  =  90°) 

(57) 


v° 

ro 


(58) 
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The  6D  expansion  coefficients  are  calculated  as, 


Voooo(rj  z) 
Vooio  (r j  z) 
Voon  (r,  z) 
V2OOO  (f,  z,  9) 
^2010  (r5  Z,  9) 
V2011  (r,z,6) 
V2eio  (r,  z ) 


yfA 

4 

yfA 

4 

yfA 

8 

4 

yfA 

4 

yfA 

~8~ 

yfA 

2 


Woo t{r,  z)  +  V00 h(r,  z)  +  2 V00b(r,  z)} , 

[Voo t(r,z)  -  V00 h(r,  z)\ , 

[Voo t(r,  z)  +  Voo h{r,  z)  -  2Vmb (r,  z)] , 

[V2ot(r,  0)  +  Vm(r,  2,  Q)  +  2V20 b(r,  z,  61)] , 
[V20 t(r,z,9)  -  V20h(r,z,9)] , 

[V20 t(r,  z,  9)  +  V20h(r,  z,  9)  -  2 V20b(r,  z,  61)] , 
V2 eb(r,  z), 


(59) 


where  A  is  the  area  of  the  surface  unit  cell,  and  the  spherical  harmonics  Ylm(9,  <p)  are 
defined  for  degree  l  =  0  and  order  m  —  0, 


yo  1  /I 

0  2  V  tt’ 


(60) 


for  degree  l  —  2, 


sin2  0  exp(2i0), 


>72(M) 


sin2  0  exp(— 2i(p), 


(61) 
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and  for  a  combination  of  F22  and  Y2  2  [174], 


Y2e(d,cf>) 


[y22(M)  +  r2-2(M)], 


15 

167T 


sin2  6  cos  20, 


where  i  is  the  imaginary  number  i  =  y  —  1. 

Finally,  the  6D  potential  Vqq  is  constructed  as, 


Vk>(z,  y ,  z,  r,  6,  0)  =H00(x,  y)  Voooofc  z)  F0° 

+  H00(x,y)V2000(r,z,6)Y2°(6) 

+  Hw(x,  y)  Vboi0(r,  z)  F0° 

+  H10(x,y)V2010(r,z,6)Y2°(6) 

+  Hu(x,y)  10011  (r,  2)  F0° 

+  Hn(x,  y)  10011  (r,  2,  0)  Y^°(0) 

+  HBll0(x,y)  V2el0{r,  z)  Y£(6,  0), 


where  the  plane- wave  diffraction  functions  are  defined  as  [86], 


H00(x,y ) 
H10(x,y) 
#n0>2/) 
HBlio(x,y) 
G 


1 

A’ 

T 


cos  Go;  +  cos  Gy] , 


2\  /  —  [cos  Gx  x  cos  Gy] , 


[cos  Gx  —  cos  Gy] , 


2tt 
az  ’ 


(62) 


(63) 


(64) 
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Note  that  in  Equation  (63),  V®  is  not  a  function  of  A. 


3.3  Classical  Trajectory  (CT)  Formulation 

The  trajectories  of  desorbing  molecules  are  modeled  classically  (CT)  with  the 
aid  of  Hamiltonian  mechanics.  The  equations  of  motion  are  derived  for  a  linear  ro¬ 
tor  molecule,  and  energies  are  identified  as  translational,  rotational,  or  vibrational. 
Quantum  effects  are  ignored,  a  necessary  assumption  for  engineering  applications. 
Otherwise,  computational  times  would  be  prohibitively  long  for  large-scale  simula¬ 
tions. 

The  kinetic  energy  T  of  a  linear  rotor  molecule  of  (total)  mass  M  and  reduced 
mass  /i  is  formed  by  taking  the  sum  of  two  components:  (1)  the  motion  of  the  center 
of  mass  about  the  origin,  and  (2)  the  motion  of  the  atoms  about  the  center  of  mass. 
In  mathematical  terms  [164:318,  330], 

T  =  ( x 2  +  y2  +  z 2)  +  ^/i  ^r2  +  r20 2  +  r202  sin2  9  j  ,  (65) 

where  ( x ,  y,  z )  is  the  location  of  the  center  of  mass  in  the  Cartesian  reference  frame, 
(r,  6 ,  0)  is  the  location  of  the  two  atoms  in  spherical  coordinates  with  respect  to  the 
center  of  mass  (r  is  the  interatomic  distance,  6  is  the  elevation,  and  0  is  the  azimuth), 
and  the  dot  notation  indicates  a  derivative  with  respect  to  time  t  (see  Section  3.6.1 
for  a  further  explanation  of  these  coordinates).  There  are  six  degrees  of  freedom 
since  the  molecule  has  two  atoms.  The  second  term  in  Equation  (65)  can  be  derived 
by  noting  that  for  a  single  particle  (mass  m)  described  by  spherical  coordinates,  its 
velocity  is  v  =  r  r  +  rO  6  +  ref)  sin  9  cf>  (r,  0 ,  and  cf>  are  the  unit- vectors  in  the  r-,  6-, 
and  0-directions,  respectively),  its  kinetic  energy  is  v  ■  vm/2,  and  the  linear  rotor 
system  is  obtained  by  replacing  m  with  /i.  Note  that  for  a  homonuclear  system  of 
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mass  m,  T  is  only  comprised  of  the  motion  of  the  center  of  mass  about  the  origin,  so 
that 

T  =  ^ m  ( x 2  +  y2  +  z 2)  .  (66) 

The  Lagrangian  L  —  T  —  V,  in  terms  of  the  generalized  coordinates  x,  y,  z,  r,  9, 
and  0,  and  the  generalized  velocities  x,  y ,  z,  r,  9,  and  0,  is  written  as 

C  =  (x2  +  y1  +  i2)  +  i/i  ^r2  +  r202  +  r202  sin2  9  j  —  V,  (67) 


where  E  is  the  potential  energy  of  the  PES.  The  generalized  momenta  are  then 
calculated  from 

dC 
Pi  =  jp~ 
oqi 

where  Pi  and  g*  are  the  generalized  momenta  and  generalized  velocities,  respectively, 
so  that 


dC 

Px  =  -WT  =  Mx, 

OX 


dC 

Pr  =  W  = 


DC  nr- 
=  My’ 


dc  H 

Pe=m=pr<>’ 


8C  nr- 
Pz  =  -wr  =  Mz, 
oz 


dC 


P4> 


=  /ir20  sin2  9 


(69) 


The  generalized  velocities  can  then  be  expressed  as  x  =  px/M,  y  =  py/M ,  z  =  pz/M, 
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r  =  pr/  p,  9  =  pe/{pr2)1  and  (j)  =  p^/^pr2  sin2  9).  Thus,  T  can  be  rewritten, 


T  = 


,y^2  I  ~»2  I  fy-.Z  ‘ 

Px+Py+Pz  ,  P 


2  M 


+^+A+ 


p i 


2/i  2/xr2  2/ir2  sin2  61 


+  Pg  + 

2  M 


v2  72 

+  —  H - 

2/i  2/  ’ 


(70) 


where  /  =  /ir2  is  the  reduced  moment  of  inertia,  and  the  rotational  angular  momen¬ 
tum  J"  is  defined  by 


r-=Pt  +  A 


sin2  9 


(71) 


The  Hamiltonian  Pi  (or  total  energy)  can  now  be  calculated  from  the  Legendre 
transformation  of  C  [131], 


7L  =  Y PiQi  -  A 

i 


=  pxx  +  Pyij  +  pzz  +  prf  +  pg9  +  p^d)  -  (T  -  V) , 


Px  +  Pi  +  Pz 

M 


Pr 

+  —  + 


P2e 


+ 


p\ 


p  pr* 


pr2  sin2  9 


T+ V, 


=  2  T  -  T  +  V, 


=  r  +  y, 


Px  +  P2y  +  P2z 

2  M 


v2  72 
+  Pjl  +  2_  +  V. 

2/i  2/ 


(72) 
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Equation  (72)  could  have  also  been  obtained  from  the  knowledge  that  since  the 
molecule  is  acting  under  a  conservative  force  (i.e.  a  closed  system),  %  =  T  +  V . 

The  rotational  angular  momentum  J  from  Equation  (71)  is 


J2 


P2e  + 


Pi 

sin2  6 ' 


(71) 


From  the  Schrodinger  equation  for  a  rigid  rotor,  we  have  that  J2 
the  rotational  number  J  is 


J  = 


-1+./Z7I 

2  V  n2  4’ 


J(J+l)fi2.  Thus 
(73) 


where  h  =  h/ (2n)  is  the  reduced  Planck  constant,  and  h  is  the  Planck  constant.  Note 
that  J  must  be  calculated  with  the  dimensional  value  of  J . 

The  vibrational  number  nv  is  calculated  from  the  vibrational  term  of  the  Hamil¬ 
tonian  Hv  =p2/(2/i)  in  Equation  (72), 


nv 


Hv_  1 

hv  2’ 


(74) 


where  v  is  the  vibrational  frequency.  The  value  of  v  for  H2  in  the  gas  phase  is  taken 
here  to  be  v  =  1.242  x  10 14  cycles  per  second  (see  Table  1  in  Reference  [131]). 

Hamilton’s  equations  then  uniquely  describe  the  time  evolution  of  the  molecule, 


Qi  = 


&H 

dpt  ’ 


(75) 
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so  that 


and 
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dU 

dV 
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(76) 


(77) 


3.4  Non-Dimensionalization 


Since  simulations  are  conducted  on  the  molecular  scale,  equations  are  non-dimen- 
sionalized  to  avoid  unnecessary  numerical  errors.  The  non-dimensional  parameters 
reported  here  are  not  available  elsewhere  in  the  literature.  Therefore,  this  formulation 
is  novel  and  would  improve  the  accuracy  of  current  CT  implementations  if  employed. 

Non-dimensionalization  is  achieved  with  the  aid  of  the  mass  of  the  molecule 
M,  the  reduced  mass  of  the  molecule  /x,  the  length  of  the  surface  unit  cell  a;,  the 
most-probable  speed  of  a  gas-phase  molecule  under  equilibrium  conditions  ump  = 
\j2kTs/M  =  yj 2RspTs ,  a  characteristic  energy  kTs,  the  most  probable  time  that 
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it  takes  one  gas-phase  molecule  under  equilibrium  conditions  to  cross  the  unit  cell 
imp  =  ai/ump,  two  different  characteristic  linear  momenta  (Mump  and  /nimp),  and 
one  characteristic  angular  momentum  aipump.  Here,  k  is  the  Boltzmann  constant, 
Ts  is  the  temperature  of  the  surface,  and  Rsp  =  k/M  is  the  specific  gas  constant.  The 
non-dimensional  parameters  (denoted  by  an  asterisk  superscript)  are 
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x  =  —, 
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V 
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ft  ^mp 
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Dmp 


(78) 


Note  that  even  though  6  and  (j)  are  already  non-dimensionalized  (radians),  they  are 
included  for  completeness. 

With  these  parameters,  the  non-dimensionalized  Hamiltonian  PL*  is 


2  *2  *2  P  *2  P  3 


n / *  * z  ,  *z  ,  *z  ,  r 

H  =Px  +Pv  +Pz  +  77  P 


M  M  r 


*2 


+  W\ 


(79) 
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where  the  non-dimensionalized  rotational  angular  momentum  J*  is 


J 


*2 


*2  ,  P* 
Pe  +~ 


sin2  9*  ’ 


(80) 


and  the  non-dimensionalized  Hamilton’s  equations  of  motion  are 


and 
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(82) 


where  the  dot  notation  indicates  a  derivative  with  respect  to  the  non-dimensionalized 
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time  t* ,  and 
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^mp 


Ump 
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et 


mpi 


and 


<fi*  =  4>tmp,  (83) 
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(84) 


For  the  case  where  the  atoms  have  equal  mass  (i.e.  m  =  mi  =  m2),  p  =  m/2  and 
p/M  =  1/4.  This  quantity,  p/M,  is  the  product  of  the  non-dimensionalized  masses 
of  each  atom, 

p  mi  m2 

M  =  ~M  ~M' 

and  as  such  could  be  considered  another  non-dimensional  parameter. 


3.5  Keck’s  Method 

To  calculate  a  trajectory,  one  must  supply  initial  conditions.  Intuitively,  one  might 
set  the  initial  conditions  of  the  adsorbate,  then  integrate  the  equations  of  motion 
forwards  in  time  until  (and  if)  the  molecule  desorbs.  However,  it  was  found  early  on 
that  by  applying  initial  conditions  in  this  manner,  only  a  small  number  of  trajectories 
would  even  reach  the  reaction  zone,  or  TS.  Keck  found  an  easier  way.  Basing  his  work 
on  Transition  State  Theory  (TST)  [89;  139],  Keck  noted  that  every  valid  reaction  must 
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pass  through  the  TS.  Hence,  he  chose  to  set  initial  conditions  at  the  TS,  rather  than 
outside  of  it.  He  then  integrated  the  equations  of  motion  forwards  and  backwards  in 
time,  observing  whether  or  not  the  full  trajectory  represented  a  chemical  reaction.  If 
not,  it  was  rejected.  He  noticed  that  most  of  the  trajectories  were  accepted.  Due  to 
its  success,  Keck’s  method  is  applied  here  [18;  72;  73]. 

A  couple  of  observations  can  be  made  at  this  point  about  Keck’s  method.  First, 
microscopic  reversibility  must  hold  in  order  to  integrate  the  equations  of  motion 
backwards  in  time.  Second,  an  assumption  must  be  made  about  the  initial  conditions. 
Fortunately,  it  has  been  shown  that  an  equilibrium  distribution  is  present  at  the  TS 
when  products  are  absent  [3].  For  desorption,  an  absence  of  products  implies  that 
desorption  is  occuring  in  a  vacuum,  a  valid  assumption  here. 

3.6  Transition  State  (TS)  Determination 

TS  locations  must  be  determined  before  Keck’s  method  can  be  implemented.  The 
calculation  of  TS  locations  continues  to  be  an  active  area  for  research  [111].  As  an 
example,  one  study  on  the  H2-Cu(lll)  system  found  that  all  desorption  events  had 
a  single  TS  location  in  common  [131],  made  possible  perhaps  by  the  simplified  PES 
model  (LEPS)  employed.  Meanwhile,  another  author  has  warned  against  assuming 
that  desorption  occurs  from  a  well-defined  TS  since  the  surface  provides  a  distribution 
of  varying  activation  barriers  [60]. 

With  those  studies  in  mind,  this  research  has  taken  the  following  approach  for 
the  6D  Wiesenekker  PES  on  H2-Cu(100),  similar  to  that  found  in  Reference  [131]. 
First,  it  is  assumed  that  hydrogen  is  adsorbed  to  the  surface  either  on  top  (T),  hollow 
(H),  or  bridge  (B)  sites.  Then,  the  20  unique  H-H  adsorbate  configurations  with  the 
adatom  separation  distance  r  within  the  bounds  ai  <  r  <  2a/  are  considered.  For  each 
configuration,  approximate  reaction  paths  are  calculated  with  their  corresponding  TS 
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locations  via  the  Chain  algorithm.  Since  each  configuration  usually  exhibits  a  poorly- 
defined  TS,  the  average  TS  location  is  assumed  to  be  the  TS  through  which  every 
reaction  path  passes  for  that  configuration.  Finally,  of  the  20  configurations,  only 
four  were  found  to  significantly  contribute  to  thermal  desorption  (see  Figure  13). 

The  determination  of  TS  locations  involves  multiple  coordinate  systems,  explained 
in  Section  3.6.1,  and  the  use  of  the  Hessian  matrix,  explained  in  Section  3.6.2.  The 
details  of  the  Chain  algorithm  are  presented  in  Section  3.6.3,  after  which  an  explana¬ 
tion  of  how  average  TS  locations  are  determined  for  the  H2-Cu(100)  system  is  given 
in  Section  3.6.4. 

3.6.1  Different  Coordinate  Systems. 

There  are  multiple  coordinate  systems  employed  in  CT  simulations.  They  are  in¬ 
troduced  here  to  familiarize  the  reader  with  the  different  nomenclature.  First,  there  is 
the  Cartesian  coordinate  system,  which  explicitly  describes  the  (x,  y,  z )  coordinates 
of  every  atom  in  the  molecule.  For  a  molecule  composed  of  two  atoms  with  indi¬ 
vidual  Cartesian  coordinates  (xi,yi,Zi)  and  (cc2, 2/2,  ^2),  the  molecule’s  location  and 
orientation  in  Cartesian  coordinates  is  z  1,  £2, 2/2,  ^2),  as  shown  in  Figure  10(a). 

Then  there  are  the  internal  coordinates  of  a  molecules.  Again,  for  a  molecule 
with  only  two  atoms,  its  location  and  orientation  can  be  described  by  the  Cartesian 
coordinates  of  its  center  of  mass  ( x ,  y,  z)  and  the  spherical  coordinates  of  the  atoms 
(r,  6,  </>)  with  respect  to  the  center  of  mass.  In  internal  coordinates,  the  molecule  is 
represented  by  (x,y,  z,r,9,(j>),  shown  in  Figure  10(b).  To  convert  from  Cartesian  to 
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2 


Figure  10.  Cartesian  and  internal  coordinate  systems  for  a  molecule  with  two  atoms. 
The  x-,  y-,  and  ^-coordinates  for  the  internal  system  are  at  the  center  of  mass  of  the 
molecule. 


internal  coordinates, 


x  = 
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and  to  convert  from  internal  to  Cartesian  coordinates, 
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where 


r 

Ax  =  —  sin  9  cos  0, 
r 

Ay  =  —  sin  6  sin  0, 
r 

A  z  =  -  cos  6. 


(88) 


It  is  also  be  necessary  to  convert  the  time  derivatives  between  the  two  coordinate 
systems,  (iq,  yi,  iy,  x2,  y2,  z2)  and  ( x,y,  z,f,0,q f>),  respectively.  To  convert  the  deriva¬ 
tives  from  Cartesian  to  internal  coordinates,  take  the  derivative  with  respect  to  time 
of  Equation  (86)  while  recalling  the  chain  rule, 
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\Xi  ~x2J 


(89) 


To  convert  the  derivatives  from  internal  to  Cartesian  systems,  take  the  derivative 
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with  respect  to  time  of  Equations  (87)-(88), 


x\  =  x  +  —r  sin  6  cos  <f>  +  -  r  9  cos  6  cos  0  —  -  r  0  sin  9  sin  0, 


iji  =  y  +  —  r  sin  9  sin  0  +  —  r  9  cos  9  sin  0  +  -  r  0  sin  9  cos  0, 


z\  —  z  +  -  f  cos  9  —  -  r  9  sin  9, 


x2  =  x  —  —  r  sin  6  cos  (j)  —  -  r  9  cos  9  cos  (j)  +  —  r  0  sin  9  sin  0, 


y2  —  y  —  -  f  sin  9  sin  0  —  -r  9  cos  9  sin  < 


—  —  r  0  sin  6  cos  0, 


zo  —  z - r  cos  9  H —  r  9  sin  9. 

2  2 


(90) 


Note  that  the  non-dimensionalized  forms  of  Equations  (86)- (90)  are  obtained  by  sim¬ 
ply  replacing  each  variable  with  its  non-dimensional  analog. 

Even  though  Equations  (86)-(90)  are  non-dimensionalized  in  a  straight-forward 
manner,  the  non-dimensionalizations  of  the  momenta  in  Cartesian  coordinates  must 
be  handled  more  carefully.  Each  atom  is  represented  with  its  own  set  of  coordinates. 
Since  each  atom  has  its  own  mass  (mi  or  m2),  the  generalized  momenta  must  be  non- 
dimensionalized  by  the  mass  of  the  atom  they  are  representing: 


pxi  =  m  iiq,  pyi  =  mii/i,  pzi  =  m^i, 

pX2  =  m2x2,  py2  =  m2y2,  pZ2  =  m2z2,  (91) 
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and 
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(92) 


Or,  in  terms  of  velocities, 


(93) 


where  for  equal  masses  rri!  =  m2  =  M/ 2,  and  ^Jrrii/M  =  yJm2/M  =  l/y/2. 

Another  way  to  describe  the  dynamics  of  the  molecule  is  in  the  mass-weighted 
Cartesian  (MWC)  coordinate  system  [160:196-208].  Each  of  the  Cartesian  coordinates 
are  multiplied  by  the  square  root  of  their  respective  atom’s  mass, 


{xly/m[,  yiy/rrn,  zly/ml,  x2^fm~2,  y2^/rn^,  z2yjm~2) , 


where  mi  and  m2  are  the  masses  of  the  two  atoms  in  the  molecule.  The  benefit  of  the 
MWC  coordinate  system  can  be  seen  by  considering  the  Euler-Lagrange  equations  in 
Cartesian  coordinates, 

±(K\  =  ac 

dt  l  ilq,  j  Dq,  ’  k  ' 

where  q3  are  the  generalized  coordinates  (aq,  yi,  Zi,  x2,  y2,  z2),  and  q3  are  the  general¬ 
ized  velocities  (iq,  yi,  iq,  x2,  y2,  z2),  and  the  Lagrangian  C  =  T  —  V  is  composed  of 
the  kinetic  energy  T  and  the  potential  energy  V.  The  kinetic  energy  T  is  a  function 
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of  the  q3  only, 


T=\Hmrfy  (95) 

i 

The  potential  energy  V  can  be  expressed  in  the  harmonic  approximation  as 

V  =  V0  +  i  ^2  QjHjMk,  (96) 

j,k 

where  the  linear  term  is  neglected,  Vo  is  the  potential  energy  V  evaluated  at  some 
initial  configuration,  and  Hjjk  is  the  (j,  k)  element  of  the  Hessian  matrix  H  (evalu¬ 
ated  at  the  same  initial  configuration) .  By  noting  that  T  —  T(q)  and  V  =  V(q), 
Equation  (94)  can  be  expressed  as 

d  ( dT\  _  dV 
dt  \  dq3  J  dq3  ’ 


mjQj  =  ^  H3.kqk,  (97) 

k 

where  the  double-dot  notation  indicates  a  second  derivative  with  respect  to  time  t. 
However,  Equation  (97)  can  be  further  simplihed  by  expressing  it  in  MWC  coordi¬ 
nates, 

=  (98) 

k 

where  /3j  is  the  jth  element  of  the  generalized  coordinates, 

{x\ y/m{,  ziy/rnl,  x2y/rn2,  2/2 z2sJrn2) , 


and  the  MWC  Hessian  IT  has  elements 


Hj,k 

yjmjmk ' 
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(99) 


For  the  case  where  all  of  the  atoms  in  the  molecule  have  the  same  mass  (i.e.  rn3  =  rrik), 
or  for  the  case  where  the  molecule  only  has  two  atoms,  H'  is  just  a  contant  multiple 
of  H.  The  MWC  generalized  momenta  p\  can  be  expressed  in  terms  of  the  momenta 
in  Cartesian  coordinates  Pi,  so  that, 


Pxi 

/  =  p to 

Pyi  yfrn? 

,  Pz  1 

Pzi  vW 

Px  2 

/  Py2 

/  Pz2 

Py2  ~  y/m? 

Pz2  y/rm' 

(100) 


The  non-dimensionalized  MWC  momenta  p'*  are  then, 


V 


/* 

Xi 


V 


/* 

X2 


P'xx 

/  * 

Pyi  = 

Pyi 

/* 

Pz  1  = 

P'ZI 

UmpVM1 

^mp\/  -M- 

M 

to 

/* 

PV2  = 

P'V2 

/* 

Pz 2  = 

^3 

to 

UmpVM1 

^mp\/  -M- 

(101) 


It  is  interesting  to  note  that  the  non-dimensional  values  of  the  Cartesian  and  the 
MWC  momenta  are  identical  (e.g.  p*Xl  =  p'x*). 


3.6.2  Hessian  Matrix. 

The  Hessian  matrix  H  is  an  important  quantity  because  it  arises  in  the  quadratic 
term  of  a  Taylor  series  expansion.  For  the  Taylor  series  expansion  of  the  potential 
energy  V  of  the  6D  PES  [160:196-208], 


V(q0  +  A  q)  =  V(q0) 


x  ^  dv 

k  dqk 


Aqk  +  9  A(ijH^  AQk  + 


(102) 


Qo 


j,k 
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where  the  Ag-terms  represent  steps  of  the  Cartesian  coordinates 


A q  =  (Axi,  A?/i,  Azi,  Ax2,  A y2,  Az2),  (103) 


qo  is  some  initial  system  configuration,  and  H^'k  is  defined  as  the  (j,  k )  element  of 
the  Hessian  H  evaluated  at  q0 , 


d2V 


dqjdqk 


QO 


(104) 


Thus,  the  Hessian  contains  information  concerning  the  second  derivatives,  and  is 
the  multi-dimensional  analog  to  the  one- dimensional  second  derivative.  The  non- 
dimensional  form  for  a  Hessian  element  in  Cartesian  coordinates  H*k  is 


Hh  = 


kT, 


Hi,j. 


(105) 


H  is  a  symmetric  matrix,  so  there  is  a  unitary  matrix  U  that  diagonalizes  H, 
(UTHU)M  =  4  ;Afc,  where  4,;  is  the  delta  function,  Xk  is  an  eigenvalue  of  H.  and 
the  T  superscript  denotes  the  matrix  transpose.  In  addition,  the  columns  of  U  are 
the  eigenvectors  of  H.  For  non-linear  molecules,  3 N  —  6  of  the  eigenvalues  will  be 
non-zero,  while  for  linear  molecules,  3 N  —  5  of  the  eigenvalues  will  be  non-zero.  For 
a  linear  molecule  with  two  atoms  ( N  =  2),  only  one  eigenvalue  will  be  non- zero,  and 
five  will  be  zero.  Those  five  trivial  eigenvalues  correspond  to  the  translational  and 
rotational  motions  of  the  molecule.  If  the  non-zero  eigenvalue  is  actually  negative, 
then  the  location  qo  +  Aq  is  a  saddle  point  on  the  PES.  Hence,  the  problem  of  finding 
a  saddle  point  on  the  PES  (which  can  indicate  a  TS)  is  equivalent  to  locating  where 
the  Hessian  has  a  negative  eigenvalue. 

To  discretize  Equation  (104),  a  double-sided  finite  difference  is  employed.  Since 
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the  Hessian  is  a  symmetric  matrix,  Hfrk  =  To  reduce  numerical  error,  each 

element  is  of  the  Hessian  is  calculated  by  averaging  the  two  equal  values,  H^’k  and 
so  that 


dV 

dqj 

dV 

q0  +  Aqh/2 

dV 

qo-Aqk/2  ,  ^Qk 

dV 

qo  +  Aqj/2  ^qk 

qo-Aqj/2 

2Aqk 

1 

2Aqj 

(106) 


where  the  configuration  q0  +  Aqq/2  is  the  initial  configuration  q0  perturbed  only  in 
the  j-direction  by  the  value  Aqj/2,  and  the  first-order  derivatives  are  calculated  with 
a  central  finite  difference  scheme, 

V  (go  +  +  A q^j  -  V  (q0  +  -  Aq^j 

qo  +  Aqj/2 


dV 

dq. 


V  I  q0  +  -  A qj  -F  go  -  -Ag 


2A  qj 


(107) 


and 


dV_ 

dqj 


V  q0 


A  qj 


A  q3 


v  q0 


A  qj 


A  q3 


qo-A.qj/2 


2A  q,j 


v  \q0  +  -  A qj  -F  go  -  xA q 


2Aqj 


(108) 


The  non-dimensionalized  form  of  Equation  (106)  is  obtained  by  simply  starring  each 
of  the  quantities. 

When  assigning  initial  values  at  the  TS,  it  is  important  to  consider  the  eigenmode 
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basis  of  H.  Mathematically,  the  eigenmode  basis  can  be  seen  by  employing  the  unitary 
matrix  U  in  Equation  (102), 

V(q0  +  A q)  =  V(q0)  +  g%. Aq  +  ^AqTU0Aq  H - , 

=  V(q0)  +  9o UUrAq  +  -2AqTXJXJTHoUUTAq  +  •  •  •  , 

=  V(Qo)  +  GqQ  +  -QrA0Q  +  •  •  •  ,  (109) 


since  UU7  =  U7  U  =  I  by  definition  (I  is  the  identity  matrix),  g0  is  the  gradient  of 
V  at  qo  with  components 


%  = 


dV 

dqk 


(110) 


QO 


Ao  is  a  diagonal  matrix  with  the  eigenvalues  of  Ho  as  elements  of  the  diagonal, 
Q  =  UTAg  is  the  component  of  the  step  Aq  along  the  eigenvector  basis  set,  and 
G0  =  Ur<7o  is  the  component  of  the  gradient  g0  along  the  eigenvector  basis  set.  For  a 
linear  rotor.  H'  is  a  constant  multiple  of  H.  and  therefore  the  unit-length  eigenvectors 
of  H'  and  H  are  identical. 


3.6.3  Chain  Algorithm. 

The  determination  of  TS  locations  and  reaction  paths  are  active  research  areas  [30; 
142],  and  fall  under  the  more  general  Mountain  Pass  Theorem  [64],  The  Chain 
algorithm,  developed  by  Liotard  and  Penot  [94;  95],  addresses  a  specific  need  to 
calculate  the  TS  location  for  a  (approximate)  reaction  path,  described  by  a  chain  of 
discrete  points. 

The  Chain  algorithm  takes  as  inputs  the  initial  (adsorbed)  and  final  (desorbed) 
configurations,  as  well  as  the  known  PES.  An  initial  path  is  chosen  that  connects  the 
initial  and  final  configurations,  which  is  then  modified  until  the  TS  is  found.  The 
TS  is  located  at  the  saddle  point  with  the  largest  energy  value,  while  a  saddle  point 
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occurs  where  a  negative  eigenvalue  of  the  Hessian  matrix  exists.  The  (approximate) 
Hessian  is  calculated  on-the-fly  via  gradient  calculations  and  the  Broyden-Fletcher- 
Goldfarb-Shanno  (BFGS)  algorithm  [136:521-523].  Note  that  the  linear  dependence 
threshold  mentioned  by  Liotard  and  Penot  [94;  135]  is  overcome  by  instead  utilizing 
BFGS. 

The  initial  reaction  path  is  created  with  two  straight  lines  connecting  known 
starting  and  ending  configurations  M\  and  M2,  respectively,  which  are  both  min¬ 
ima  on  the  PES.  The  intermediate  point  connecting  M\  with  M2  is  taken  as  the 
starting  configuration  but  with  the  intermolecular  distance  of  the  ending  configu¬ 
ration  M2.  In  other  words,  if  the  starting  configuration  M\  in  internal  coordinates  is 
(xi,  yi,  z i,  r i,  6*i ,  0i),  and  the  ending  configuration  M2  is  ( x 3,  2/3,  Z3,  r 3,  #3,  03),  then 
the  intermediate  configuration  is  (xi,  y  1,  z1:  r3,  fa)  since  r3  is  the  intermolecular 
distance  of  M2.  For  each  line,  a  linear  interpolation  between  each  of  the  coordi¬ 
nates  creates  the  chain  of  points.  In  this  work,  each  line  was  divided  equally  into  10 
segments. 

The  calculation  steps  of  the  Chain  algorithm  are: 

1.  Set  the  threshold  parameters  ( l\  =  0.01  Z2  and  /o  =  0.01 1\). 

2.  Determine  the  point  pH  which  currently  approximates  the  TS.  This  is  done  by 
evaluating  the  potential  V  at  every  point  on  the  chain,  pn  is  the  point  with 
the  highest  value  of  V. 

3.  Calculate  the  negative  gradient  of  the  potential  at  pn,  G  =  —  W(p/f),  in 
Cartesian  coordinates.  Each  of  the  partial  derivatives  is  calculated  with  a  cen¬ 
tral  difference, 

dV  _  V(xi  +  Ax)  -  V(xi  -  Ax) 

dxi  ~  2Ax  1  j 

where  Ax  =  Iq. 
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4.  Update  the  Hessian  H  via  the  BFGS  algorithm.  If  the  point  pn  has  not  moved 
significantly  in  the  most  recent  step,  then  H  is  not  updated.  In  other  words,  if 
llPifnew  -Pf/oidll  <  tlien  Hnew  =  Hold.  For  this  work,  e,  =  10^5az. 

5.  Calculate  the  quadratic  direction,  Q  =  H  1 G  and  determine  if  Q  is  operative. 

(a)  Q  is  operative  if  the  index2  of  H  is  unity,  ||Q||  <  l2,  and  Q  ■  G  > 
e2  ||Q||  ||G||,  where  e2  =  0.1. 

6.  If  the  magnitude  of  the  negative  gradient  ||G||  falls  below  a  given  threshold 
gc  (i.e.  1 1 G 1 1  <  gc),  then  stop;  pn  is  the  TS.  The  gradient  threshold  gc,  below 
which  the  PES  is  assumed  to  be  flat,  is  set  at  0.04  kTs  az_1. 

7.  Select  the  step  direction  D  from  either  the  quadratic  direction  Q ,  the  direction 
of  the  projected  gradient  G*,  or  one  of  the  unit  vectors  A  of  the  links  adjacent 
to  pH. 

(a)  The  quadratic  direction  Q  is  chosen  if: 

i.  Q  is  operative,  Q  is  decreasing  (Q  G  >  0),  and  ||G||  <  g i,  where  g\ 
is  a  pre-convergence  threshold  on  the  gradient  norm  (gi  =  10 gc  here); 
or 

ii.  Q  is  operative,  Q  is  decreasing,  and  9  >  g,  where  9  is  the  angle 
between  G  and  the  pseudo-tangent  to  the  path  T  (9  e  [0,  90°]),  and  g 
is  a  threshold  value  ( g  =  30°  here).  The  pseudo-tangent  T  is  defined  as 
the  unit  vector  parallel  to  the  straight  line  between  Ph-i  and  Ph+i, 
or 

T=  Ph+i-Ph-1  (112) 

_  IIph+1  -  Vh-i\\ 

2The  index  of  H  is  defined  as  the  number  of  negative  eigenvalues  of  H. 
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(b)  The  projected  gradient  G*  is  defined  as  G*  =  G  —  (G-T)T.  The  direction 
of  G*  is  chosen  if: 

i.  Q  is  operative,  Q  is  increasing,  and  9  >  77;  or 

ii.  Q  is  not  operative,  Q  is  decreasing,  and  9  >  77;  or 

iii.  Q  is  not  operative  or  ||G||  <  g  1,  and  9  <g. 

(c)  The  more  increasing  unit  vector  A  of  the  two  links  adjacent  to  pn  is  chosen 

if: 

i.  Q  is  operative,  Q  is  increasing,  and  ||G||  <  gp, 

ii.  Q  is  not  operative,  Q  is  increasing,  and  9  >  77;  or 

iii.  If  neither  unit  vector  A  is  increasing,  then  no  search  direction  is  chosen 
for  this  step.  Skip  to  Step  10,  with  p*  chosen  as  the  midpoint  between 
Ph-i  and  pH+i,  or 


*  Ph-i  +  Ph+i 

= - 5 - 


(113) 


8.  Conduct  a  line  search  to  find  the  optimal  step  length  A  G  [Zo,  h\-  In  the  direction 
D ,  evaluate  V  at  10  equally-spaced  points  beginning  at  Z0  and  ending  at  h-  The 
point  with  the  lowest  value  of  V  determines  the  step  length. 

9.  Calculate  the  new  point  p *,  where 

P*  =  Ph  +  Xjn=jTr  (114) 


10.  Replace  pH  with  p* . 

11.  If  necessary,  reconstruct  the  chain  by  eliminating  meanders  and/or  inserting 
new  points. 
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12.  Decrease  the  value  of  l2.  Here,  l2nev,  =  0.9Z2okr  The  value  of  l2  for  the  first  step 

is  2  A. 


13.  Return  to  Step  1. 


The  BFGS  algorithm  (see  Step  4)  is  a  quasi-Newton  method,  since  it  utilizes  not 
the  actual  Hessian  but  an  approximation  to  it.  BFGS  assumes  that  the  PES  can  be 
approximated  (locally)  by  the  Taylor  series  expansion  in  Equation  (109),  so  that  at 
step  i, 

V(ql0  +  q)  «  V (<?o)  +  9oTq  +  ^qTH'0q.  (115) 

At  each  subsequent  step  i  +  1,  the  Hessian  Hq  is  updated  to  Hg+1  with  information 
gained  from  the  step  taken  (g(,+1  —  ql0),  and  from  the  gradients  at  steps  i  and  i  +  1, 
which  are  gl0  and  g'a+ 1 .  respectively.  In  mathematical  terms, 


Hq+1 


0 0o  -  Q'o)  ®  (0o  -  q o)  [HJ,  ■  (gf1  -  gi,)}  ®  [HJ,  •  (gf1  -  g'0)] 


(q o  -  qlo )  •  too  -  0o ) 


(^+i  -  gi)  ■  Hq  •  0 0 o+i  -  0&) 


[(0o+  -  0o )  '  Ho  '  (0o+  -  0o)]  u®u, 


(116) 


where  the  vector  u  is  defined  as 


=  (0o+1  -  00 ) _ Ho  •  (0o+1  -  00 )  /117N 

(0o+1  -  00 )  •  (0o+1  -  0o)  (0o+1  -  0o)  •  Ho  •  (0o+1  -  0o)  ’ 

and  <g)  is  the  outer  product  of  two  vectors,  resulting  in  a  matrix.  For  example,  the 
(i,j)  component  of  v  <S>  w  is  ViWj.  The  Hessian  at  the  first  step  is  taken  to  be  the 
identity  matrix. 

To  avoid  severe  kinks  in  the  reaction  path,  the  chain  is  reconstructed  and/or  new 
points  are  inserted  in  Step  11,  as  shown  in  Figure  11.  For  the  insertion  of  new  points, 
if  either  of  the  two  links  connected  to  Ph  have  lengths  longer  than  l2,  then  a  new  point 
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;  pi+ 1  Ph 


(a)  Before 


(b)  After 


Figure  11.  Two-dimensional  representation  of  a  meander  bypass  in  the  Chain  algorithm. 
Here,  Pi+i  is  inside  the  hypersphere,  so  the  link  is  made  between  Pi+i  and  p/j.  However, 
if  pi  happened  to  be  the  point  in  the  hypersphere,  then  the  new  link  would  be  made 
between  pi  and  pjj. 

is  either  placed  halfway  through  that  link,  or  at  a  distance  Z2  from  pn,  whichever 
is  smaller.  To  bypass  meanders  in  the  chain,  a  hypersphere  is  constructed  around 
point  Ph  of  radius  rHs  (set  to  2  A).  Beginning  at  the  start  of  the  chain  (with  Mi  as 
the  first  point  of  the  first  link),  each  link  is  tested  in  order  until  pn  is  reached.  Let 
point  pi  be  the  start  of  the  link  under  consideration,  and  let  Pi+i  be  the  end  of  the 
same  link.  If  the  link  crosses  the  hypersphere,  then  the  chain  is  modified  in  one  of 
two  ways.  Either  pn  is  connected  with  pi  or  with  Pi+i-,  depending  on  whether  Pi+i 
is  outside  of  the  hypersphere  or  not,  respectively.  As  a  result,  all  of  the  points  that 
existed  on  the  chain  between  pi  (or  Pi+i)  and  pn  are  permanently  removed.  The 
same  procedure  is  repeated,  except  that  it  is  begun  at  the  end  of  the  chain  (with  M2 
as  the  last  point  of  the  last  link). 

An  example  reaction  path  calculated  with  the  Chain  algorithm  is  presented  in 
Figure  12.  The  adsorbed  configuration  corresponds  to  r  =  ai,  while  the  desorbed 
configuration  is  located  at  z  —  1.8  a/.  The  x-,  y -,  9-,  and  0-values  are  constant  in  this 
example.  The  atoms  start  in  the  adsorbed  configuration,  travel  along  the  2D  PES  to 
the  TS,  and  then  continue  on  to  the  desorbed  configuration.  The  calculated  reaction 
path  is  approximate;  output  of  the  Chain  algorithm  is  only  the  TS  location. 

The  example  shown  in  Figure  12  is  only  for  two  dimensions.  Even  though  Chain 
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Figure  12.  An  example  path  calculation  with  the  Chain  algorithm.  The  important 
data  point  is  the  TS  location  (circle).  The  approximate  reaction  path  (dashed  line)  is 
only  useful  in  arriving  at  the  TS  location,  and  is  discarded  after  the  calculations  are 
completed,  x,  y ,  9 ,  and  <j)  are  constant  in  this  example.  The  adsorbed  configuration 
is  at  r  =  ai,  while  the  desorbed  configuration  is  at  z  =  1.8  a*.  Contour  lines  are 
expressed  in  units  of  eV.  This  contour  plot  is  a  2D  slice  of  the  6D  PES  for  x  —  ai, 
y  —  ai/2,  9  —  90°,  and  <fi  —  90°,  copresponding  to  the  case  wherein  hydrogen  atoms  are 
adsorbed  on  top  of  adjacent  copper  atoms,  and  desorbed  as  H2  above  a  bridge  site. 


calculations  must  consider  all  six  dimensions,  each  link  can  be  uniquely  described  by 
five  2D  lines, 

aj  =  m]a\  +  bj,  (118) 

where  oq  is  one  of  the  six  Cartesian  coordinates  (not  neccesarily  x\ ),  j  is  an  index  that 
cycles  through  the  remaining  five  coordinates,  bj  is  the  aq-intercept  (when  aq  =  0), 
and  rrij  is  the  slope  of  the  2D  line, 


O7  2+1  ®  j,  i 

rrij  =  — - — , 

04,2+1  —  al,i 


(119) 


where  the  subscripts  i  +  1  and  i  signify  that  the  coordinates  are  from  points  p^+i  and 
Pi ,  respectively.  a±  must  be  chosen  such  that  rrij  is  defined.  However,  unless  pi  and 
Pi_ |_i  are  the  same  point,  there  will  be  at  least  one  coordinate  that  can  be  chosen  as 


The  equations  to  determine  whether  or  not  a  link  crosses  the  hypersphere  were 
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not  found  in  the  literature.  Thus,  they  were  derived  for  this  work.  The  equation  for 
the  hypersphere,  centered  at  (aqc,  a2c>  a3 c>  a4c>  ot5c,  a6c),  is 

6 

rHS  =  _  «ic)2, 

l=i 

6 

=  («!  —  «ic)2  +  —  aJC)2, 

J=2 

6 

=  (ai  -  aic)2  +  ^  [(mjCti  +  bj)  -  ajc ]2  .  (120) 

1=2 

Solving  for  ap,  the  generalized  quadratic  formula  is  obtained, 


on 


—b  ±  a/62  —  4ac 
2  a 


(121) 


where 


6  6 

6  =  rrijajC  —  ai c, 

1=2  1=2 


c 


1 

2 


6 


.1=2 


)2  + 


HS 


(122) 


Once  cti  is  calculated,  the  rest  of  the  crossing-point  coordinates  can  be  determined 
from  Equation  (118).  The  crossing  points  are  then  tested  to  determine  if  they  lie  on 
the  link  in  question  or  not. 


3.6.4  TS  Locations  for  H2-Cu(100). 

This  section  describes  in  detail  how  adsorbed  and  desorbed  configurations  are 
chosen  (inputs  to  the  Chain  algorithm),  and  how  the  most  important  (average)  TS 
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locations  are  selected  for  CT  simulations. 


The  Chain  algorithm  requires  adsorbed  and  desorbed  configurations  as  inputs. 
Adsorbate  configurations  are  assumed  to  be  a  subset  of  the  symmetry  adsorption 
sites  on  the  copper  surface.  A  symmetry  site  is  one  in  which  the  hydrogen  atom  is 
adsorbed  either  above  a  bridge  site  (B),  a  hollow  site  (H),  or  directly  on  top  of  a 
copper  atom  (T).  Since  there  are  many  unique  symmetry  configurations  on  a  surface, 
only  a  subset  of  symmetry  configurations  may  be  chosen  to  determine  the  most  likely 
TS  locations  that  contribute  to  H2  desorption.  If  only  symmetry  configurations  where 
the  two  H  atoms  are  no  less  than  a;  (but  no  more  than  2  a;)  apart,  a;  <  r  <  2 a;,  are 
considered,  then  there  are  20  unique  symmetry  configurations.  All  of  the  adsorbed 
configurations  are  initially  assumed  to  be  parallel  to  the  surface  (0  =  90°)  while 
the  height  above  the  surface  z  is  optimized  by  minimizing  V  while  keeping  x,  y,  r, 
6,  and  (j)  constant.  The  configurations  are  then  relaxed  by  minimizing  V  within  a 
hypersphere  of  radius  rns  =  0.01  a;  by  taking  an  optimal  step  in  the  direction  of  the 
negative  gradient  of  the  potential  G  =  —W  [see  Equation  (111)].  This  optimization 
strategy  is  also  referred  to  as  a  line  search.  A  line  search  is  therefore  conducted  to 
find  the  optimal  step  length  A  G  [uhs/100,  rHs)  in  the  direction  of  G.  To  do  so,  V  is 
evaluated  at  10  equally-spaced  points  beginning  at  rHs/100  and  ending  at  rHs-  If  it 
is  found  that  A  =  rns,  then  rns  is  increased  by  0.01  a;,  and  the  line  search  is  repeated 
until  r'Hs/100  <  A  <  rns-  This  technique  is  similar  to  the  trust  radius  method  [35]. 
Examples  of  initial  and  relaxed  adsorbed  configurations  are  shown  in  Figure  13. 

Desorbed  configurations  have  as  their  height  z  =  9ao  =  1.87  a;  (above  which 
the  influence  of  the  PES  is  negligible),  and  as  their  intermolecular  bond  length 
r  =  0.291a;,  which  is  the  equilibrium  bond  distance  for  an  H2  molecule  in  the  gas 
phase  [175].  The  remaining  locations  and  angles  are  chosen  randomly  with  a  uniform 
distribution  over  the  unit  cell  and  the  unit  sphere,  so  that  x  G  [0,  a;],  y  G  [0,  a/] , 
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Figure  13.  The  symmetry  sites  found  to  contribute  the  most  to  thermal  desorption  un¬ 
der  the  CT  construct,  as  seen  perpendicular  to  the  £cy-plane.  Empty  circles  represent 
the  copper  surface  atoms,  while  filled  circles  represent  adsorbed  hydrogen  atoms.  The 
numbers  correspond  to  those  found  in  Table  1.  The  perfectly  symmetric  (original)  con¬ 
figurations  are  shown  in  13(a),  where  the  adsorbed  hydrogen  atoms  are  located  either 
on  top  (T),  bridge  (B),  or  hollow  (H)  sites.  After  conducting  a  trust  radius  calculation 
on  the  original  configurations,  the  actual  (relaxed)  configurations  are  obtained  [shown 
in  13(b)].  It  is  interesting  to  note  that  the  relaxed  configurations  tend  to  position  the 
midpoints  near  bridge  sites. 
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6  e  [0,27r],  and  0  G  [0, 7r] .  Note  that  for  a  random  location  on  the  unit  sphere,  6 
can  be  assigned  a  uniform  random  angle  by  6  =  2nRi ,  while  0  must  be  assigned  a 
random  angle  through  the  relation  0  =  cos_1(2i?2  —  1),  where  R\  and  R2  are  uniform 
random  numbers  on  [0, 1]. 

Each  adsorbed  configuration  is  matched  with  5000  random  desorbed  configura¬ 
tions,  from  which  5000  TS  locations  are  determined.  By  taking  an  average  TS  location 
for  each  symmetry  configuration,  20  average  TS  locations  are  identified  as  possible 
contributors  to  thermal  desorption.  To  determine  which  TS  locations  significantly 
contribute  to  thermal  desorption,  trajectories  were  calculated  for  19  of  the  20  config¬ 
urations3  at  a  surface  temperature  of  Ts  =  1100  K.  The  desorption  angle  results  were 
then  compared  with  experimental  values  (see  Figure  16(a)),  and  each  configuration 
was  weighted  to  optimize  the  fit  between  simulation  and  experiment.  The  correspond¬ 
ing  relaxed  adsorption  configurations  are  shown  in  Table  1,  the  weights  and  average 
TS  locations  are  shown  in  Table  2,  and  a  comparison  between  TS  locations  and  re¬ 
laxed  adsorption  configurations  is  shown  in  Figure  14.  It  was  found  in  this  research 
that  only  four  configurations  were  required  to  match  desorption  angle  experimental 
values  for  that  specific  temperature.  It  is  assumed  that  those  four  configurations 
dominate  thermal  desorption  over  all  temperature  regimes. 

Trajectories  were  calculated  by  integrating  Equations  (82)  forwards  and  back¬ 
wards  in  time  from  the  TS  accroding  to  Keck’s  method.  The  ode45  function  in 
MATLAB™  solves  nonstiff  differential  equations  with  fourth-order  accuracy,  and  was 
used  to  perform  the  calculations.  It  employs  the  explicit  Runge-Kutta  (4,  5)  formula 
(the  Dormand-Prince  pair).  As  such,  ode45  is  a  one-step  solver,  meaning  that  it  only 
requires  information  from  the  immediately-preceding  step.  For  integration  forward  in 

3One  of  the  configurations  was  rejected  due  to  unreasonable  simulation  times  (more  than  four 
hours  per  trajectory  on  average).  That  configuration  corresponded  to  the  H  atoms  both  adsorbed 
at  bridge  sites  (2  a/  apart),  desorbing  as  H2  above  the  bridge  site  halfway  between  them. 
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Figure  14.  A  comparison  of  the  average  TS  locations  (black)  and  the  relaxed  adsorption 
configurations  (grey),  as  seen  perpendicular  to  the  ccy-plane.  The  TS  locations  tend  to 
have  midpoints  directly  above  the  midpoints  for  the  relaxed  adsorption  configurations. 

Table  1.  Relaxed  adsorption  symmetry  site  values. 


Site  Label 

Internal  Coordinates 

x  [ ai ] 

y  M 

2  [ai] 

r  [ai] 

e 

0 

1 

0.50 

0.00 

0.59 

2.00 

O 

o 

o 

0° 

2 

0.05 

0.55 

0.52 

1.57 

o 

o 

C} 

00 

o 

3 

0.00 

0.60 

0.48 

1.84 

o 

o 

O} 

170° 

4 

0.45 

0.00 

0.56 

1.50 

o 

o 

o 

0° 

time  from  the  TS,  integration  stops  when  the  molecule’s  height  above  the  surface  is 
z  =  1.87  cii ,  which  is  the  criterion  for  desorption.  If,  however,  the  molecule  dissociates 
(r  >  0.784  ai )  before  desorbing,  then  the  trajectory  is  rejected  [131].  For  integration 
backward  in  time  from  the  TS,  integration  stops  once  the  molecule  has  dissociated. 
Similarly,  the  molecule  cannot  desorb  before  it  dissociates  in  the  backwards  direction 
in  order  for  the  trajectory  to  be  considered  valid. 


83 


Table  2.  Average  TS  locations  with  their  respective  weights. 


TS  Label 

Weight  [%] 

Internal  Coordinates 

x  [ai\ 

y  M 

z  k] 

r  [ai] 

e 

0 

1 

44.6 

0.42 

0.04 

0.71 

0.50 

96° 

176° 

2 

32.2 

0.01 

0.53 

0.52 

0.45 

115° 

28° 

3 

12.2 

0.99 

0.60 

0.50 

0.35 

o 

o 

00 

o 

4 

11.0 

0.46 

0.01 

0.60 

0.54 

100° 

165° 

3.7  Initial  Conditions 

With  Keck’s  method,  a  trajectory  begins  at  the  TS  and  is  calculated  forwards 
and  backwards  in  time.  Therefore,  initial  conditions  correspond  to  the  momenta  at 
the  TS  location.  The  momenta  are  found  by  first  selecting  a  random  value  of  the 
molecule’s  energy  from  a  truncated  MB  distribution  (adjusted  for  the  energy  lost  by 
the  adsorbate  in  overcoming  the  activation  barrier),  and  then  randomly  distributing 
that  energy  among  the  six  momenta  in  the  MWC  Hessian  eigenspace  with  equal 
probability.  The  MWC  momenta  in  the  Hessian  eigenspace  are  then  converted  to 
MWC  momenta  in  the  Cartesian  eigenspace,  after  which  the  momenta  are  converted 
to  internal  coordinates  for  calculation  purposes. 

3.7.1  Truncated  Maxwell-Boltzmann  (MB)  Distribution. 

The  adatoms  are  assumed  to  be  in  equilibrium  at  the  surface  temperature,  de¬ 
scribed  by  the  MB  distribution  fE, 

hdE=2 V?fe)  exp  (“ w) iE’  (123) 
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where  /#  dE  is  the  probability  that  an  adatom  will  have  an  energy  in  the  range  E  to 
E  +  dE.  In  non-dimensional  units, 


f  E*  dE*  =  2 


exp  (~E*)dE*. 


(124) 


However,  once  the  molecules  arrive  at  the  TS,  they  have  lost  some  energy.  This  energy 
is  simply  the  difference  in  potential  energy  between  the  adsorbed  configuration  and 
the  TS  configuration,  AH.  The  energy  distribution  of  the  molecules  at  the  TS  can 
then  be  described  by  a  truncated  MB  distribution  /Tg, 


(125) 


where  Equation  (124)  has  been  shifted  to  the  left  by  AV*  =  A V/(kTs),  and  only 
positive  values  of  E*  are  valid.  Discrete  values  are  chosen  randomly  from  this  distri¬ 
bution  in  an  accept-reject  manner  [20:423-428],  which  requires  the  distribution  to  be 
normalized  by  its  maximum  value.  Equation  (124)  has  its  maximum  at  E*  =  1/2, 
while  Equation  (125)  peaks  at  E*  =  1/2  —  AV*  with  a  value  of  /max, 


(126) 


However,  this  maximum  only  occurs  for  AH*  <  1/2.  For  AH*  >  1/2,  the  maximum 
occurs  at  E*  =  0  with  the  value 


(127) 
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Therefore,  the  fully  normalized  and  truncated  MB  distribution  /^s,  shown  in  Fig¬ 
ure  15,  is 


( 


E*  +  AH*  -  i 


[2(E‘  +  AV'*)]1/2exp 


(128) 


/  E*  \  1 /2 
(aW  +  1)  exp(-£-)- 


Ay*  > 

2 


Following  the  accept- reject  method,  given  AV*,  a  random  value  of  E*  is  chosen  from 
a  uniform  distribution  E*  G  [0,  1,  where  E^ ax  is  a  discrete  approximation  for 

positive  infinity.  For  E^ax  =  12,  the  cumulative  distribution  function  (CDF)  value 
at  is  at  most  only  0.000025  less  than  unity,  so  =  12  is  the  chosen  upper 
limit.  A  uniform  random  number  is  then  selected,  R  G  [0, 1].  If  R  <  /^s( E *), 
then  E*  is  accepted.  Otherwise,  E*  is  rejected,  and  another  value  of  E*  is  randomly 
selected  until  a  value  is  accepted.  Note  that  in  the  limit  as  AV*  — >  oo,  the  normalized 
distribution  approaches  an  exponential  decay  curve,  f^s  — >  exp (-E*).  On  a  practical 
level,  for  AV*  >  10,  the  normalized  distribution  can  reasonably  be  approximated  as 


/ts  =  exp  (-E*). 


3.7.2  Mass- Weighted  Cartesian  (MWC)  Hessian  Eigenspace. 

With  E*  chosen,  the  initial  momenta  in  the  MWC  Hessian  eigenspace  p'*{  are 
determined  by  dividing  the  energy  among  the  (squared)  six  momenta  randomly  from 
uniform  distributions.  The  subscript  u  indicates  that  the  momenta  are  in  the  direc¬ 
tions  of  the  MWC  Hessian  eigenvectors  Ui.  Since  the  change  in  the  potential  energy 
has  already  been  taken  into  account,  E*  is  equivalent  to  the  non-dimensional  kinetic 


Figure  15.  Representations  of  the  fully-normalized  and  truncated  MB  distribution 
for  different  values  of  AF*  [see  Equation  (128)].  The  shaded  area  is  the  actual 
distribution,  starting  at  E*  —  0.  The  dashed  portion  of  the  curve  shows  where  the  MB 
distribution  has  been  truncated,  and  has  as  its  minimum  abscissal  value  -AF*. 


energy  of  the  molecule  T*,  which  can  be  written  as 


t'  =  Y.pj^ 


i=  1 


For  six  random  uniform  numbers  E*  6  [0, 1]  such  that, 


^>*  =  1, 


i=l 


the  momenta  can  be  expressed  as, 


Y.pZ2  =  Y.R-T'- 


i= 1 


i= 1 


such  that, 


P 


=  ±VW’ 


(129) 


(130) 


(131) 


(132) 
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If  Uj  is  the  Hessian  eigenvector  with  the  negative  eigenvalue,  then  the  corresponding 
momentum  p'*k  should  be  in  the  positive  direction, 

Ki  =  (133) 

since  that  is  the  direction  of  recombination  across  the  TS.  The  other  five  momenta, 
however,  retain  the  ±  symbol  in  Equation  (132). 

The  momenta  in  MWC  coordinates  are  then  recovered  by  multiplying  the 
momenta  in  the  MWC  Hessian  eigenspace  p'*^  by  the  Hessian  eigenvectors  Uj , 

6 

(i34) 

3= 1 

where  Uj(i)  is  the  i-th  element  of  the  eigenvector  Uj.  Momenta  in  internal  coordinates 
can  then  be  obtained  by  the  transformations  discussed  previously  [see  Equations  (76), 
(78),  (90)-(91),  and  (100)-(101)]. 

3.8  Dynamic  Simulation  of  Desorption  Events 

As  time  passes  in  a  real  gas-surface  system,  the  coverage  on  the  surface  changes  due 
to  adsorption,  desorption,  and  permeation.  Desorption  modeling  has  been  addressed 
for  each  event.  However,  one  must  calculate  when  each  desorption  event  occurs.  Also, 
the  rate  of  atoms  arriving  at  the  surface  via  permeation  must  be  identified.  Therefore, 
Section  3.8.1  introduces  the  topic  of  permeation  and  how  it  can  be  applied  to  the 
current  model,  and  Section  3.8.2  applies  MC  event  timing  to  the  DSMC  framework. 

3.8.1  Permeation. 

As  hydrogen  molecules  desorb  from  the  surface,  new  hydrogen  atoms  are  supplied 
to  the  surface  via  permeation  at  the  rate  Q.  When  a  pressure  differential  exists, 


molecules  tend  to  diffuse  from  areas  of  high  concentration  to  areas  of  low  concen¬ 
tration,  even  in  the  presence  of  a  membrane  or  structural  component,  such  as  a 
magnetron  wall.  For  the  example  of  a  magnetron,  the  outer  casing  experiences  atmo¬ 
spheric  pressure,  while  the  inner  cavity  is  held  at  a  sub-atmospheric  pressure.  Thus, 
molecules  from  the  atmosphere  tend  to  permeate  the  structure  because  of  the  pressure 
differential,  and  find  their  way  to  the  surface  of  the  inner  cavity.  Hydrogen  perme¬ 
ates  easily  through  many  materials,  including  copper.  The  process  of  permeation 
(and  then  desorption)  for  hydrogen  on  copper  is  [54]: 

1.  The  (atmospheric)  hydrogen  gas  first  adsorbs  to  the  outer  surface,  dissociating 
into  hydrogen  adatoms  in  the  process. 

2.  Hydrogen  adatoms  are  then  absorbed  into  the  copper  bulk  according  to  the 
hydrogen-copper  solubility. 

3.  Hydrogen  atoms  then  diffuse  through  the  copper  bulk  to  the  inner  surface, 
driven  by  the  pressure  differential  existing  between  the  two  surfaces,  and  deter¬ 
mined  by  the  hydrogen-copper  diffusivity. 

4.  Upon  reaching  the  inner  surface,  hydrogen  atoms  become  adatoms,  adsorbed 
to  the  surface. 

5.  Through  a  desorption  process,  the  hydrogen  molecule  desorbs  from  the  copper 
surface  into  the  low-pressure  environment. 

By  inspection,  it  is  apparent  that  Steps  1  and  2  are  the  reverse  processes  of  Steps  4  and  5 
Permeation  is  a  function  of  the  pressures  at  both  surfaces,  the  solubility  and  diffusiv¬ 
ity  of  the  system,  as  well  as  the  temperature  of  the  copper.  These  dependencies  will 
become  more  apparent  as  the  relevant  mathematical  relations  are  presented. 


The  total  permeation  Q  of  molecules  per  unit  time  t  through  the  membrane  is 
calculated  by 

Q  =  J  ■  A,  (135) 

where  J  is  the  flux  of  molecules  traveling  through  the  area  A,  and  the  dot  product 
is  necessary  so  that  only  the  flux  normal  to  the  surface  is  considered. 

Fick’s  First  and  Second  Laws  describe  how  the  flux  J  and  concentration  of 
molecules  0  varies  within  the  material  according  to  diffusion,  in  steady-state  and 
general  conditions,  respectively  [54;  153].  Fick’s  Second  Law,  analogous  to  the  heat 
equation,  is 

^  I )  V2o.  (136) 

where  0  is  the  concentration,  D  is  the  diffusion  coefficient,  t  is  time,  and  V2  is 
the  Laplace  operator.  Under  steady-state  conditions,  Fick’s  Second  Law  reduces  to 
Laplace’s  equation, 

V20  =  0,  (137) 

which  can  be  further  manipulated  to  obtain  Fick’s  First  Law, 

J  =  —D  V0,  (138) 


where  V  is  the  gradient  operator. 

Boundary  conditions  at  both  surfaces  are  required  to  solve  Equation  (138),  which 
are  determined  via  Sievert’s  Law.  This  law  states  that  the  concentration  of  atoms 
0  dissolved  in  a  metal  is  directly  proportional  to  the  gas-solid  solubility  S  and  the 
square  root  of  the  fugacity  /  of  the  gas  molecules  in  thermodynamic  equilibrium, 

<t>  =  Sy/f,  (139) 
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where  the  fugacity  /  is  equal  to  the  gas  partial  pressure  P  for  an  ideal  gas.  With 
the  assumption  that  the  desorption  surface  is  under  vacuum  conditions  (P  =  0), 
Equation  (138)  can  be  solved  for  a  given  geometry.  For  a  planar  geometry, 


A<S>\/P 


(140) 


where  d  is  the  thickness  of  the  membrane,  P  is  the  non- vacuum  pressure,  and  4>  =  DS 
is  the  permeability,  while  for  a  cylindrical  geometry, 


2irh  <h  y/P 


(141) 


Q 


where  h  is  the  height  of  the  cylinder,  a  is  the  outside  radius,  and  b  is  the  inside 
radius.  Because  of  the  assumption  of  vacuum  conditions,  Equations  (140)-(141)  give 
the  maximum  permeation  rate.  Real  conditions  will  therefore  exhibit  a  lower  Q-value. 

The  permeability  $  is  identified  with  the  product  of  diffusivity  D  and  solubility 
S,  which  is  independent  of  the  specific  geometry, 


$  =  DS. 


(142) 


Typically,  D  and  S  are  expressed  in  terms  of  temperature  by 


(143) 
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where  Dq  and  So  are  constants,  Ed  and  E$  are  activation  energies,  k  is  the  Boltzmann 
constant,  and  T  is  temperature.  Therefore,  the  permeability  $  is  expressed  as 

*(r)  =  *0expf-^y  (144) 

where  <h0  =  D0S0  and  E$  =  ED  +  Eg-  Values  of  the  constants  and  activation  energies 
for  Hydrogen  in  various  materials  can  be  found  in  Reference  [153]. 

3.8.2  Timing. 

The  timing  of  events  is  an  important  part  of  simulating  dynamic  desorption  and 
adsorption  [46;  61;  92],  The  DSMC  algorithm  simulates  the  discrete  Bernoulli  process. 
By  associating  a  time  step  with  each  event,  DSMC  desorption  produces  a  Markov 
chain  of  events  [61],  which  simulates  a  Poisson  process  since  the  system  is  large  and 
contains  only  independent  events  [46]. 

The  probability  density  of  times  between  successive  events  ft  for  a  Poisson  process 
is  [46] 

ft  =  rexp(-rt),  (145) 

with  the  average  time  between  events  (t)  being 

(t)  =  1  (146) 

where  r  is  the  event  rate.  Instead  of  invoking  the  average  time  between  events  as  the 
time  step,  however,  time  is  incremented  after  each  desorption  event  by 

At=--lnR,  (147) 

r 

where  R  is  a  uniform  random  number  between  0  and  1,  and  r  may  be  different  for 
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each  event.  Equation  (147)  differs  from  the  corresponding  equation  in  Reference  [46] 
by  determining  At  only  from  the  desorption  rate  of  that  particular  event,  rather  than 
from  the  sum  of  the  rates  of  all  events  from  all  subpopulations.  The  distinction  is 
made  here  since  time  is  incremented  in  a  different  fashion  than  it  is  with  MC  methods. 
In  MC,  the  time  step  is  incremented  after  a  certain  number  of  events  have  passed 
(e.g.  50),  while  in  these  simulations,  time  is  incremented  after  every  event. 

The  total  rate  of  desorption  r  is  calculated  from  the  Polanyi-Wigner  model  [77], 

r  =  f  =  -■'"*“»>(-#)  <148> 

where  v  is  a  frequency  factor,  x  is  the  reaction  order  for  desorption,  E  is  the  aver¬ 
age  desorption  activation  energy,  and  N  is  the  number  of  adatoms  on  the  surface. 
Values  of  u,  x,  and  E  are  published  in  the  literature.  Anger  et  ah  present  these 
values  for  H2-Cu(lll)  and  H2-Cu(110)  [5].  However,  since  the  surface  of  Cu(100) 
undergoes  a  reconstruction  at  high  hydrogen  coverage  [31],  H2-Cu(100)  TPD  spectra 
are  published,  but  values  for  v,  x,  and  E  are  not.  Fortunately,  their  TPD  spectra 
can  be  analyzed  independently  by  the  Threshold  TPD  (TTPD)  method  presented  in 
References  [57;  77;  117]. 

To  calculate  r  for  the  next  time  step  {i  +  1),  values  must  be  taken  from  step  ( i ) 
(e.g.  Q,  r,  and  At),  since  r  is  a  function  of  the  number  of  adatoms  found  on  the  surface, 
supplied  by  permeation  and  adsorption,  and  depleted  via  desorption.  Therefore,  the 
number  of  adatoms  available  for  desorption  at  time  step  (i  +  1)  is  calculated  as 

iV(m)  =  [Q«(tW)  -  r(i)]  A t(i),  (149) 

where  the  superscript  indicates  the  time  step  from  which  values  are  taken.  Substi¬ 
tuting  Equation  (149)  into  Equation  (148),  r  for  the  (i  +  l)th  time  step  is  readily 
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calculated  by 


(150) 


r(i+1)  =  -v  [lV(i+1)]xexp 


E 

kT^+l) 


During  a  DSMC  simulation,  each  desorption  event  should  be  determined  in  the 
following  manner.  Before  time  step  (i  +  1)  begins,  if  the  time  that  has  passed  since  the 
last  desorption  event  is  greater  than  A t<yl\  then  desorption  should  occur  at  the  local 
surface  temperature  T^.  A  new  value  of  At  should  be  calculated  at  each  time  step 
until  this  condition  is  reached,  unless  if  the  surface  temperature  remains  unchanged 
from  the  previous  time  step. 


3.9  Results 

Results  from  the  current  method  are  presented  in  this  section.  Emphasis  is  placed 
on  comparing  theoretical  simulations  with  the  most  up-to-date  experimental  data 
available  for  thermal  desorption  from  the  H2-Cu(100)  gas-surface  system  [158].  These 
results  are  for  a  surface  temperature  of  Ts  =  1030  K  and  include  state-resolved  mean 
translational  energy  values,  state-resolved  translational  energy  distributions,  and  the 
quadrupole  alignment  parameter  Aq  ’ .  In  addition,  the  desorption  angle  distribution 
at  Ts  =  1100  K  is  shown.  Recall  from  Section  3.6.4  that  contributing  TS  locations  and 
their  respective  weights  were  determined  by  fitting  simulation  results  to  experimental 
data. 


3.9.1  Desorption  Angle  Distribution. 

In  order  to  compare  theory  with  experimental  results  from  Reference  [12],  a  geo¬ 
metrical  factor  based  on  the  experimental  measuring  equipment  must  be  taken  into 
account.  Simulation  results  give  a  desorption  angle  for  each  molecule,  while  exper¬ 
imental  data  provide  the  number  of  molecules  that  arrive  at  a  detector  covering  a 
finite  solid  angle  D,  at  multiple  locations  in  space. 
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As  an  example,  for  the  measurements  reported  in  Reference  [157],  the  detector 
is  9  mm  away  from  the  surface  with  a  diameter  of  2  mm,  giving  a  solid  angle  of 
=  0.0384  steradians,  whereas  an  entire  sphere  covers  47t  steradians.  Thus,  the 
detector  covers  only  a  portion  (a  cap)  of  the  surface  of  the  sphere,  while  all  simulated 
molecules  with  identical  desorption  angles  have  the  same  latitude.  To  compare  the 
two,  a  latitudinal  band  is  constructed  by  rotating  the  cap  about  the  surface  normal. 
This  band  has  a  finite  width  corresponding  to  the  diameter  of  the  detector.  The  cap 
covers  an  area  Acap  of 

Acap  =  27ri?2(l  —  cos  Ad),  (151) 

where  R  is  the  distance  of  the  detector  from  the  surface,  Ad  =  tan-1  (r/R),  and  r 
is  the  radius  of  the  detector.  To  continue  this  example,  R  =  9  mm,  r  =  1  mm, 
A 6  =  6.34°,  and  Acap  =  3.11  mm2.  The  area  of  the  latitudinal  band  Aband  is  found 
to  be 

Aband  =  47t R2  sin  9  sin  A 9,  (152) 

where  6  is  the  desorption  angle  of  interest,  measured  from  the  surface  normal.  How¬ 
ever,  Equation  (152)  only  holds  when  6  >  A 9.  For  the  case  where  9  <  A 9,  Equa¬ 
tion  (152)  must  be  modified  to  account  for  overlapping  areas, 

Aband  =  27 tR2  (1  +  sin  9  sin  A 9  —  cos  9  cos  A 9) .  (153) 


To  then  compare  simulation  results  with  experimental  data,  one  needs  to  only  multi¬ 
ply  the  simulation  values  by  the  ratio  Acap/Aband, 


A 


cap 


A 


band 


'  1  —  cos  A  9 

1  +  sin  9  sin  A 9  —  cos  9  cos  A 9  ’ 
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9  >  A 9. 


(154) 
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Figure  16.  Distribution  of  desorption  angle  0  with  respect  to  the  surface  normal  for 
H2-Cu(100)  at  Ts  =  1100  K.  Simulation  results  are  compared  with  experiment  and 
a  curve  fit  to  experimental  data.  The  sampling  is  from  500k  trajectories.  Dashed 
curve,  cosd(0)  with  d  =  5;  □,  experimental  data;  — ,  simulation  results.  Experimental 
data  are  from  Reference  [12],  Figure  3.  Differs  from  Figure  17  by  (Acap/Tband)  — 1  in 
Equation  (154). 


With  the  aid  of  Equation  (154),  desorption  angle  results  are  compared  in  Fig¬ 
ure  16(a)  to  experimental  values  taken  at  Ts  =  1100  K.  As  discussed  in  Section  3.6.4, 
contributing  TS  locations  were  determined  by  optimizing  the  fit  of  simulation  results 
at  Ts  =  1100  K  to  the  experimental  data.  Figure  16(a)  illustrates  how  well  desorption 
data  can  be  duplicated  by  this  method  and  the  Wiesenekker  PES.  The  dashed  curve 
represents  a  curve  fit  to  experimental  data  of  the  form  cos d{6)  with  d  =  5.  A  cosine 
distribution  (d  =  1)  would  occur  with  complete  thermalization  of  the  molecules  to  the 
surface  temperature  if  there  were  no  surface  corrugation.  However,  a  corrugated  PES 
alters  the  distribution  from  that  of  a  pure  cosine.  Thus,  even  though  the  molecule  is 
assumed  to  be  fully  thermalized  to  the  surface  temperature,  the  a  cosine  distribution 
is  still  not  observed  due  to  the  PES. 

Simulation  results  match  well  the  experimental  data  points  in  Figure  16(a),  ex¬ 
cept  for  desorption  angles  below  10.5°.  A  desorption  peak  is  predicted  at  2°  from 
the  normal,  a  feature  not  shown  by  the  experimental  measurements.  Further  des- 
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Table  3.  Filtering  ratios  for  desorption  angles  0  <  10.5°. 


Range  [°] 

Filtering  Ratio 

Range  [°] 

Filtering  Ratio 

[0.5, 1.5) 

0.773 

[5.5,  6.5) 

0.644 

[1.5,  2.5) 

0.456 

[6.5,  7.5) 

0.802 

[2.5,  3.5) 

0.467 

[7.5,  8.5) 

0.798 

[3.5, 4.5) 

0.515 

[8.5,  9.5) 

0.942 

[4.5,  5.5) 

0.594 

[9.5,10.5) 

0.938 

orption  experiments  would  need  to  be  conducted  to  determine  whether  or  not  that 
peak  is  physical.  There  is  no  indication  in  the  literature  that  such  a  peak  has  been 
found  before.  However,  measurements  are  typically  reported  at  large  intervals  due 
to  limitations  in  the  experimental  equipment  (i.e.  10°,  20°,  etc.),  so  that  if  a  peak 
does  exist  near  2°,  it  would  not  have  been  found  unintentionally.  However,  the  peak 
appears  to  be  non-physical  and  simply  represents  a  limitation  of  the  current  method. 
One  could  mitigate  this  peak  in  engineering  applications  by  randomly  rejecting  excess 
molecules.  To  illustrate  this  possibility,  consider  a  molecule  desorbing  with  an  angle 
between  1.5°  and  2.5°.  Since  the  experimental  curve-fit  value  at  6  =  2°  is  45.6%  of 
the  predicted  value  in  Figure  16(a),  randomly  accept  (reject)  the  molecule.  Choose  a 
uniform  random  number  R  between  0  and  1.  If  R  <  0.456,  then  accept  the  trajectory. 
Otherwise,  reject  it  and  calculate  a  new  trajectory.  The  desorption  angle  ranges  and 
ratios  of  experimental  to  predicted  values  are  listed  in  Table  3.  Note  that  since  the 
results  in  Figure  16(a)  are  binned  with  a  bin  width  of  1°,  the  desorption  angle  ranges 
have  a  width  of  1°  and  are  centered  at  integer  values. 

Figure  16(b)  shows  the  modified  desorption  angle  results  with  the  accept-reject 
criteria  taken  into  account  for  desorption  angles  6  <  10.5°.  Less  than  4%  of  the 
trajectories  are  rejected.  Even  though  it  appears  that  a  much  higher  percentage  of 
molecules  would  need  to  be  filtered  when  comparing  Figures  16(a)  and  16(b),  recall 
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Figure  17.  Raw  distribution  of  desorption  angle  6  with  respect  to  the  surface  normal 
for  H2-Cu(100)  at  Ts  =  1100  K.  Simulation  results  are  compared  with  experiment  and  a 
curve  fit  to  experimental  data.  The  sampling  is  from  500k  trajectories  (unfiltered)  and 
481k  (filtered).  Dashed  curve,  cosd(0)  with  d  —  5;  □,  experimental  data;  — ,  simulation 
results.  Experimental  data  are  from  Reference  [12],  Figure  3.  Differs  from  Figure  16 
by  Acap/Aband  in  Equation  (154). 


that  simulation  values  are  modified  per  Equation  (154)  before  they  are  compared  with 

experimental  data.  Thus,  the  predicted  raw  count  is  only  modified  by  less  than  4%. 

This  idea  is  illustrated  in  Figure  17,  where  simulation  results  are  presented  in  their  raw 

count,  while  experimental  results  are  modified  by  in  Equation  (154). 

Figures  16  and  17  contain  identical  data,  just  presented  in  a  different  manner.  The 

filtering  effect  in  Figure  17(b)  is  only  for  desorption  angles  6  <  10.5°,  where  it  is  clear 

that  only  a  small  portion  (less  than  4%)  of  the  trajectories  are  filtered  out. 

The  remaining  figures  will  be  presented  with  filtering  taken  into  account.  Energy 

distributions  see  less  than  a  2%  change  (see  Figures  18  and  19)  and  less  than  a 

9%  change  (see  Figure  20)  when  filtering  is  applied.  More  variability  is  seen  in 

the  translational  energy  distributions  than  in  the  mean  translational  energy  values 

(2  \ 

due  to  mathematical  averaging.  Likewise,  the  quadrupolc  alignment  parameter  Tig 
experiences  less  than  a  3%  change  over  its  entire  range  (see  Figures  21,  22,  and  23). 
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3.9.2  Energy  Distributions. 


The  average  translational  energy  of  desorbing  H2  is  shown  in  Figure  18  as  a  func¬ 
tion  of  the  rotational  number  J  for  two  different  vibrational  numbers  (nv  =  0,1), 
where  J  and  nv  are  calculated  by  rounding  to  the  nearest  integer.  Simulations  pre¬ 
dict  lower  average  kinetic  energies  than  is  observed  in  experiments  (0.28  eV  for  nv  =  0, 
and  0.15  eV  for  nv  —  1).  This  result  is  to  be  expected  due  to  how  the  calculations  are 
conducted.  Recall  that  the  Wiesenekker  PES  was  constructed  from  MD  simulations 
of  hydrogen  already  adsorbed  to  the  copper  surface.  However,  in  the  experiments  the 
hydrogen  adlayer  was  populated  via  permeation  through  the  bulk.  Thus,  either  the 
PES  needs  to  be  constructed  taking  into  account  the  absorption  energy  barrier,  or 
vibrational-phonon  coupling  should  be  modeled  to  more  accurately  predict  desorption 
energies.  Figure  18(a)  presents  the  simulation  data  without  the  absorption  barrier 
energy,  while  Figure  18(b)  includes  the  absorption  barrier  energy.  Predictions  match 
experiment  with  excellent  agreement.  The  one  questionable  data  point  is  at  J  =  3, 
where  experimental  measurements  show  a  sudden  drop  in  the  mean  translational  en¬ 
ergy.  Sementa  et  al.  [158]  provide  no  explanation  for  the  abnormal  data  point,  so 
that  one  is  left  to  conclude  that  it  was  an  unexpected  measurement. 

Since  simulations  are  able  to  recreate  the  shape  of  the  experimental  curves  in 
Figure  18(b)  for  both  vibrational  and  rotational  numbers,  the  absorption  barrier  is 
shown  here  to  primarily  affect  the  translational  energy  of  the  desorbing  gas,  and  then 
only  as  a  function  of  the  vibrational  energy.  For  engineering  applications,  one  can 
take  an  upper-limit  approach,  assuming  that  the  absorption  barrier  contribution  is 
a  constant  value  for  all  vibrational  numbers.  Since  there  is  no  vibrational-phonon 
coupling  for  nv  =  0,  the  value  found  here  (0.28  eV)  can  serve  as  that  upper  limit. 
The  absorption  barrier  is  equivalent  to  the  solubility  energy  E$  in  Equation  (143), 
and  the  value  of  0.28  eV  falls  within  the  published  range  for  E$  [153].  Vibrational- 
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phonon  coupling  may  explain  why  the  absorption  barrier  energy  contribution  is  less 
for  nv  —  1  than  for  nv  =  0. 

Since  the  absorption  barrier  energy  mainly  affects  the  mean  translational  energy, 
it  is  reasonable  to  conclude  that  adatoms  receive  this  additional  energy  in  a  similar 
manner.  On  average,  the  adatoms  receive  the  additional  energy  in  the  x-,  y -,  and 
^-directions  in  equal  proportion  to  each  other.  In  other  words,  the  absorption  barrier 
energy  accelerates  both  adatoms  in  the  same  direction  before  they  arrive  at  the  TS 
location. 

Some  researchers  employ  GW,  which  rejects  trajectories  that  have  vibrational 
energies  far  from  quantum  values,  in  order  to  achieve  more  realistic  statistics  from 
their  Monte  Carlo  simulations  [23;  24;  131].  However,  it  was  found  in  this  research 
that  GW  is  not  required. 

Figure  18(b)  is  a  key  result  of  this  work.  One  of  the  original  goals  of  modeling 
thermal  desorption  in  DSMC  was  to  accurately  predict  internal  and  translational 
energy  distributions.  Up  until  now,  no  research  effort  has  been  able  to  accomplish  this 
feat,  especially  without  GW.  Therefore,  this  work  provides  an  engineering  alternative 
to  time-intensive  MD  and  DFT  calculations. 

The  average  total  energy  is  presented  in  Figure  19,  shown  as  a  function  of  both 
nv  and  J.  The  average  total  energy  of  desorbing  molecules  actually  increases  with 
increasing  J.  This  result  indicates  that  two  desorbing  molecules  with  the  same  initial 
total  energy  will  transfer  different  energy  amounts  to  the  surface,  depending  on  their 
rotational  energy.  Thus,  molecules  with  lower  (final)  rotational  energies  tend  to 
transfer  more  energy  to  the  surface  than  those  with  higher  (final)  rotational  energies. 

The  translational  energy  distributions  for  ( nv  =  0,  J  =  1  —  9)  are  shown  in  Fig¬ 
ure  20.  As  the  rotational  number  increases,  the  distribution  changes  from  unimodal 
(J  =  1  —  2),  to  bimodal  (J  =  3  —  7),  and  then  back  to  unimodal  (J  =  8  —  9).  In 
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(a)  Without  Barrier  Energy  (b)  With  Barrier  Energy 

Figure  18.  State-resolved  mean  translational  energy  of  desorbing  H2  from  Cu(100)  at 
Ts  —  1030  K.  The  sampling  is  from  a  total  of  1M  trajectories  (600k  for  nv  —  0,  213k  for 
nv  =  1).  Simulation  results  are  shown  with  filled  symbols  (•  for  nv  —  0,  ▲  for  nv  —  1), 
while  experimental  values  are  shown  in  open  symbols  (o  for  nv  —  0,  A  for  nv  —  1),  taken 
from  Reference  [158],  Figure  9.  Simulations  predict  a  lower  mean  translational  energy 
due  to  the  energy  barrier  a  hydrogen  atom  experiences  when  migrating  from  the  bulk 
to  the  surface  of  copper.  The  simulation  curves  are  therefore  adjusted  upwards  in  (b) 
to  illustrate  this  effect  (0.28  eV  for  nv  —  0,  and  0.15  eV  for  nv  —  1).  Vibrational- 
phonon  coupling  explains  why  the  energy  barrier  contribution  is  less  for  nv  —  1  than 
for  nv  =  0. 
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Figure  19.  State-resolved  total  energy  of  desorbing  H2  from  Cu(100)  at  Ts  —  1030  K 
from  simulations.  The  sampling  is  from  600k  trajectories  for  nv  =  0  (•),  and  213k 
trajectories  for  nv  —  1  (A).  Total  number  of  trajectories  for  all  vibrational  numbers  is 
1M.  Total  energy  increases  for  increasing  J,  but  decreases  for  increasing  nv. 


comparing  Figures  18  and  20,  the  average  translational  energy  of  the  distributions  is 
seen  to  decrease  as  J  increases  for  most  of  the  rotational  numbers.  The  only  excep¬ 
tion  is  from  J  =  8  to  J  =  9,  where  the  average  translational  energy  increases  with 
increasing  J. 


3.9.3  Rotational  Alignment. 

The  predicted  average  rotational  alignment  of  the  molecule  as  it  desorbs  is  also 

compared  with  experiment.  The  average  rotational  quadrupole  alignment  parameter 
(2) 

Aq  ;  is  calculated  in  reference  to  the  surface  normal.  For  a  given  rotational  number 
J,  Aq2)  is  defined  as  [39;  131;  182:Chap.  5,  Application  13], 

Ao\j)  =  (3cos2S  -  l)j,  (155) 
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Figure  20.  Translational  energy  distributions  for  states  (nv  =  0,  J  —  1  —  8)  of  H2 
desorbing  from  Cu(100)  at  Ts  —  1030  K  from  simulations.  Taken  from  a  sampling 
of  1M  trajectories.  The  rotational  number  J  and  the  number  of  trajectories  at  each 
state  are  listed  in  the  subcaptions.  The  absorption  barrier  energy  is  not  included  in 
simulations.  As  the  rotational  number  J  increases,  the  distributions  transition  from 
unimodal  (J  =  1  —  2)  to  bimodal  (J  —  3  —  7)  and  then  back  to  unimodal  (J  =  8  —  9). 
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where  5  is  the  angle  between  the  £-axis  and  the  rotational  angular  momentum  vector 


J, 


cos 


Jz 

J’ 


(156) 


Jz  is  the  ^-component  of  the  angular  momentum,  and  the  brackets  (•)  j  denote  taking 

(2) 

the  average  for  all  molecules  that  have  that  value  for  J  once  desorbed.  /Lq  ;  is  valid  in 
the  range  A ^  e  [—1,  2],  with  —1  indicating  that  the  molecules  desorb  in  a  completely 


cartwheel  motion,  2  indicating  that  they  desorb  with  a  helicopter  motion,  and  any 
other  value  indicating  a  combination  of  these  two  motions. 

The  ^-component  Jz  is  easily  calculated  in  Cartesian  coordinates.  The  angular 
momentum  for  a  molecule  is  the  sum  of  the  cross-products  of  the  radius  vector  r%  and 
the  momentum  vector  pt  for  each  atom  i, 


n 


(157) 


where  n  is  the  number  of  atoms  in  a  molecule  (2  in  this  case),  and  the  angular 
momentum  for  the  ith  atom  is 


(158) 


Ji  =  rtx  pi , 


where  ry  begins  at  the  molecule’s  center  of  mass  and  terminates  at  atom  i,  and  pi  is 
the  linear  momentum  of  the  ith  atom.  The  ^-component  of  J  is  then 


n 


(159) 


where  the  ^-component  of  angular  momentum  for  the  ith  atom  Jz  i  is 


(160) 
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Figure  21.  Quadrupole  alignment  parameter  Aq  '  as  a  function  of  the  translational 
energy  of  the  desorbing  gas  H2  from  Cu(100)  for  states  ( nv  —  0,  J  —  4)  and  ( nv  = 
1,  J  =  2  —  4)  at  Ts  —  1030  K.  The  sampling  is  from  1M  total  trajectories.  The  solid 
lines  represent  simulation  results,  while  the  dashed  lines  represent  experimental  values 
from  Reference  [158],  Figure  6.  The  lines  are  both  calculated  with  a  ten-point  moving 

average  of  the  raw  data.  In  (b),  the  average  value  of  Aq  ;  is  given  over  the  rotational 

•  •  •  •  •  •  (2)  • 
numbers  considered.  Helicopter  motion  is  indicated  by  A],  ’  —  2,  and  cartwheel  motion 

(2) 

by  Aq  =  —  1.  The  horizontal  solid  line  is  for  reference  only.  Simulation  values  match 
experiment  fairly  well. 


where  rXi  =  Xi  —  x  and  ryi  =  Hi  —  y  [see  Equations  (87)-(88)],  and  pXi  and  pyi  are  the 
x-  and  ^-components  of  momentum  for  the  ?'th  atom,  respectively. 

The  quadrupole  alignment  parameter  Aq  ;  provides  a  measure  of  how  well  the 
physics  are  represented  by  the  simulation.  It  indicates  how  the  molecules  are  oriented 
as  they  rotate  in  space.  Even  though  Aq  ;  may  not  be  beneficial  in  engineering  codes, 
it  is  discussed  here  to  validate  the  model.  As  shown  in  Figure  21,  the  model  performs 
well  when  compared  against  experiment.  The  predicted  curves  follow  experiment 
qualitatively,  and  are  within  15%  and  12%  for  {nv  =  0,  J  =  4)  and  (n„  =  1,  J  = 
2  to  4),  respectively.  The  molecules  tend  to  desorb  with  more  of  a  cartwheel  than 
helicopter  motion. 

Figure  22  shows  the  average  quadrupole  alignment  parameter  Aq  for  the  ground 
vibrational  state  nv  =  0  at  various  rotational  numbers.  Except  for  J  =  8,  predictions 
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(b)  J  =  2,  32k  (c)  J  =  3,  33k 


Translational  Energy  [eV]  Translational  Energy  [eV] 

(d)  J  =  5,  32k  (e)  J  =  8,  63k 


(2) 

Figure  22.  Quadrupole  alignment  parameter  Aq  ;  as  a  function  of  the  translational 
energy  of  the  desorbing  gas  H2  from  Cu(100)  for  states  ( nv  —  0,J  =  1,2,  3,  5,  8)  at 
Ts  —  1030  K.  The  sampling  is  from  a  total  of  1M  trajectories.  The  rotational  number  J 
and  the  number  of  trajectories  are  given  in  each  subcaption.  The  solid  lines  represent 
simulation  results,  while  the  dashed  lines  represent  experimental  values  from  Refer¬ 
ence  [158],  Figure  14.  The  lines  are  both  calculated  with  a  ten-point  moving  average 

of  the  raw  data.  Helicopter  motion  is  indicated  by  A g2^  =  2,  and  cartwheel  motion  by 
(2) 

Ag  =  —  1.  The  horizontal  solid  line  is  for  reference  only.  Simulation  values  match 
experiment  very  well,  except  for  a  portion  of  the  plot  from  state  ( nv  —  0,  J  —  8). 
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(a)  J  =  2,  7k  Trajectories  (b)  J  =  3,  10k  Trajectories 


(c)  J  =  4,  10k  Trajectories 

Figure  23.  Quadrupole  alignment  parameter  A],  as  a  function  of  the  translational  en¬ 
ergy  of  the  desorbing  gas  H2  from  Cu(100)  for  states  ( nv  —  1,  J  —  2  —  4)  at  Ts  —  1030  K. 
The  sampling  is  from  a  total  of  1M  trajectories.  The  rotational  number  J  and  the  num¬ 
ber  of  trajectories  are  given  in  each  subcaption.  The  solid  lines  represent  simulation 
results,  while  the  dashed  lines  represent  experimental  values  from  Reference  [158],  Fig¬ 
ure  15.  The  lines  are  both  calculated  with  a  ten-point  moving  average  of  the  raw  data. 
Helicopter  motion  is  indicated  by  Aq  ’  —  2,  and  cartwheel  motion  by  A],  '  =  —  1.  The 
horizontal  solid  line  is  for  reference  only.  There  is  a  decent  match  between  simulation 
values  and  experiment. 
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agree  well  with  experimental  values.  The  maximum  difference  between  the  two  are 
7%  (J  =  1),  8%  (J  =  2),  10%  (J  =  3),  7%  (J  =  5),  and  22%  (J  =  8).  Cartwheel 
motion  is  favored  for  the  most  part.  For  J  =  8,  the  model  predicts  more  pronounced 
cartwheel  motion  than  experiment. 

In  general,  the  trend  appears  to  favor  cartwheel  motion  for  low  translational  energy 
values.  As  the  translational  energy  increases,  the  motion  on  average  has  an  equal 
contribution  of  cartwheel  and  helicopter  motion.  An  explanation  for  this  may  be  that 
at  low  translational  energies,  the  desorbing  molecules  are  influenced  much  more  by 
the  PES  than  at  high  translational  energies.  As  the  molecule  desorbs,  it  experiences  a 
longer  residence  time  on  the  surface  for  a  low  translational  energy.  Therefore,  the  PES 
has  a  greater  influence  on  the  molecule’s  orientation.  On  the  other  hand,  a  molecule 
with  a  high  translational  energy  spends  a  shorter  amount  of  time  being  influenced 
by  the  PES,  and  so  its  orientation  is  much  more  dependent  on  the  molecule’s  initial 
conditions. 

For  nv  —  1,  Figure  23  presents  Ay0  as  a  function  of  translational  energy.  This 
vibrational  state  is  less  populated  than  the  ground  state.  Hence  the  statistics  in 
Figure  23  have  a  higher  variability  than  is  seen  in  Figure  21.  In  spite  of  the  fact  that 
fewer  trajectories  are  represented  in  Figure  23,  general  trends  can  still  be  assessed. 
Predictions  follow  experimental  data  well  qualitatively,  even  though  quantitatively 
there  is  some  deviation.  The  model  shows  a  maximum  difference  from  experiment  of 
12%  for  J  =  2,  19%  for  J  =  3,  and  14%  for  J  =  4. 

3.10  Chemisorption  Summary 

A  new  and  successful  model  for  thermal  desorption  was  developed  in  this  chap¬ 
ter.  It  can  be  directly  incorporated  into  Direct  Simulation  Monte  Carlo  codes  as  a 
boundary  condition.  Appropriate  timing  has  been  introduced  as  well  so  that  non- 
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equilibrium,  temporally-accurate  simulations  can  be  conducted.  A  typical  activated 
chemisorption  system,  H2-Cu(100),  was  used  to  demonstrate  the  novel  method.  How¬ 
ever,  this  work  is  not  limited  to  hydrogen  nor  to  copper.  It  is  expected  that  thermal 
desorption  due  to  the  Langmuir-Hiushelwood  reaction  on  any  gas-surface  system  can 
be  modeled  with  the  same  techniques. 

At  least  seven  significant  contributions  to  the  state-of-the-art  have  been  devel¬ 
oped  in  this  research.  Thermal  desorption  modeling  is  now  available  for  engineering 
applications,  to  include  desorption  angle,  internal  energies,  translational  energy,  and 
the  molecular  alignment.  Accurate  temporal  modeling  in  Direct  Simulation  Monte 
Carlo  has  been  introduced.  The  equations  of  motion  have  been  presented  in  a  non- 
dimensionalization  scheme,  something  not  found  in  the  literature.  The  new  model  is 
not  limited  by  the  number  of  transition  states  it  can  consider.  In  fact,  this  research  has 
shown  that  many  transition  states  should  be  identified  to  provide  accurate  results. 
Initial  conditions  were  determined  from  a  truncated  Maxwell- Boltzmann  distribu¬ 
tion,  developed  here  with  its  accompanying  accept-reject  form.  Gaussian  weighting, 
a  filtering  scheme  which  greatly  increases  run  times,  has  been  shown  here  to  be  un¬ 
necessary  under  the  proposed  model.  Finally,  this  research  has  clearly  shown  that 
the  absorption  barrier  energy  from  permeation  not  only  significantly  contributes  to 
the  translational  energy  of  desorbing  molecules,  but  it  also  has  little  effect  on  their 
rotational  and  vibrational  energies. 
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IV.  Conclusions 


Surface  processes  are  varied  and  complex,  and  exhibit  unique  behaviors  for  differ¬ 
ent  gas-surface  systems.  There  is  a  need  to  be  able  to  simulate  these  processes  for  a 
wide  variety  of  systems  with  engineering  methods.  This  work  has  sought  to  address 
that  need  by  applying  the  Direct  Simulation  Monte  Carlo  framework  to  the  areas  of 
adsorption  and  desorption  in  physisorptive  and  chemisorptive  systems,  respectively, 
something  that  has  not  been  accomplished  until  now. 

For  adsorption,  the  Xe-Pt(lll)  system  was  chosen  as  representative  of  a  typical 
physisorptive  system.  Building  on  the  Cercignani-Lampis-Lord  scattering  kernel  and 
the  Multi-Stage  scattering  model,  a  new  method,  the  Modified  Kisliuk  with  Scattering 
method,  was  successfully  developed  for  calculating  adsorption  probabilities  and  scat¬ 
tering  properties.  The  results  show  that,  for  the  Xe-Pt(lll)  system,  this  method  was 
able  to  accurately  predict  adsorption  probabilities  as  a  function  of  coverage,  even  for 
a  gas-surface  system  where  adsorption  probabilities  increase  with  increased  coverage. 

The  Modified  Kisliuk  with  Scattering  method  is  applicable  to  any  gas-surface 
physisorption  system.  Thus,  even  though  the  Xe-Pt(lll)  system  was  investigated 
here,  any  combination  of  a  gas  and  surface  can  be  applied,  as  long  as  the  bond  is 
physisorptive  in  nature. 

One  could  also  extend  this  method  to  include  non-physisorptive  systems  by  as¬ 
suming  that  no  chemical  reactions  occur.  In  other  words,  one  could  assume  that 
gas  molecules  only  scatter  or  adsorb  in  relation  to  a  surface,  without  the  presence  of 
dissociation,  recombination,  or  other  chemical  reactions.  In  engineering  applications 
where  quick  calculations  are  required,  this  assumption  could  provide  reasonable  re¬ 
sults.  Or,  a  comparison  with  experimental  results  could  illustrate  the  effects  due  to 
chemical  reactions  alone.  In  this  manner,  the  contributions  made  by  surface  chemical 
effects  could  be  uncoupled  from  those  made  by  scattering  and  adsorption. 
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The  Modified  Kisliuk  with  Scattering  method  dramatically  improves  the  scattered 
translational  energy  profile.  Experimental  data  show  that  for  some  scattering  angles, 
the  incident  molecules  actually  gain  energy.  The  Cercignani-Lampis-Lord  kernel  alone 
cannot  predict  this  behavior.  However,  with  the  new  method,  this  super-elastic  scat¬ 
tering  behavior  can  be  predicted. 

The  scattering  angle  distribution  is  also  improved  by  the  Modified  Kisliuk  with 
Scattering  method  over  the  traditional  Cercignani-Lampis-Lord  kernel.  Experimental 
distributions  of  Xe-Pt(lll)  are  tighter  than  the  kernel  can  predict.  The  new  method 
easily  tightens  the  scattering  distribution,  providing  more  realistic  predictions. 

For  the  first  time,  the  Cercignani-Lampis-Lord  kernel  was  applied  here  to  scat¬ 
tering  off  of  an  adlayer.  This  work  has  shown  that  precursor-mediated  gas-adlayer 
interactions  can  be  predicted  by  the  kernel,  as  incorporated  into  the  new  method. 
Previously,  work  with  the  Cercignani-Lampis-Lord  kernel  has  been  confined  to  the 
gas  scattering  off  of  the  primary  surface. 

Optimal  values  for  the  accommodation  coefficients  were  found  with  the  new  method 
by  comparing  the  predictions  with  experiment.  The  Modified  Kisliuk  with  Scattering 
method  thus  provides  an  additional  mean  for  determining  these  coefficients,  for  use 
both  with  the  Cercignani-Lampis-Lord  kernel,  as  well  as  other  scattering  models. 

For  thermal  desorption,  the  H2-Cu(100)  system  was  chosen  as  the  prototypical 
chemisorptive  system  with  activated  desorption.  Classical  trajectory  methods  on  a 
six-dimensional  potential  energy  surface  were  incorporated  with  Keck’s  method  to 
efficiently  calculate  desorption  parameters  (desorption  angle,  molecular  alignment, 
and  translational,  rotational,  and  vibrational  energies).  The  equations  were  non- 
dimensionalized  in  order  to  reduce  numerical  error.  Relevant  transition  states,  found 
with  the  Chain  algorithm,  were  weighted  according  to  their  relative  contributions  to 
the  desorption  angle.  Initial  conditions  at  the  transition  states  were  determined  with 
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a  truncated  Maxwell-Boltzmann  distribution.  A  method  for  conducting  temporally- 
accurate  simulations  is  developed,  by  taking  into  account  appropriate  event  timing 
and  permeation. 

Many  contributions  to  the  state-of-the-art  were  developed  in  this  work.  First, 
non-dimensionalized  classical  trajectory  equations  were  derived  here.  Nowhere  in  the 
literature  were  the  equations  of  motion  found  to  be  non-dimensionalized.  This  fact 
is  surprising  since  non-dimensionalization  can  greatly  reduce  numerical  error. 

Previous  researchers  have  limited  Keck’s  method  to  only  one  transition  state. 
However,  this  work  has  shown  that  it  is  beneficial  to  consider  all  contributing  transi¬ 
tion  states.  Each  location  contributes  in  a  different  manner  to  the  final  results,  with 
unique  angular  and  energy  distributions.  This  novel  approach  proved  here  to  improve 
classical  trajectory  modeling. 

The  trunctated  Maxwell-Boltzmann  distribution  was  developed  here,  along  with 
its  accept-reject  form.  The  final  trajectories  are  heavily  dependent  upon  initial  con¬ 
ditions.  Therefore,  the  truncated  Maxwell-Boltzmann  distribution  is  necessary  for 
Keck’s  method. 

Temporally-accurate  desorption  modeling  also  requires  correct  event  timing.  In 
the  literature,  it  has  been  known  for  many  years  how  to  correctly  time  desorption 
events  in  Monte  Carlo  simulations.  However,  Direct  Simulation  Monte  Carlo  has 
never  before  incorporated  this  technique.  With  a  few  modifications,  event  timing  has 
been  developed  as  a  part  of  this  research  for  use  with  Direct  Simulation  Monte  Carlo. 

It  is  common  in  classical  trajectory  modeling  to  require  Gaussian  weighting,  which 
weights  results  based  on  how  closely  their  classical  vibrational  energies  align  with 
quantum  values.  Unfortunately,  in  Direct  Simulation  Monte  Carlo  applications,  Gaus¬ 
sian  weighting  would  dramatically  increase  simulation  run  times.  This  work  has 
demonstrated  that  Gaussian  weighting  is  not  always  necessary  in  classical  trajectory 
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simulations.  Perhaps  Gaussian  weighting  is  not  required  here  because  of  the  other 
significant  contributions.  In  other  words,  it  is  possible  that  the  physics  are  so  well 
modeled  in  this  research,  that  one  does  not  need  to  take  a  quasi-classical  trajectory 
approach  here. 

Finally,  there  has  been  some  discussion  in  the  literature  as  to  how  much,  or  in 
what  manner,  the  absorption  barrier  energy  contributes  to  thermal  desorption.  This 
debate  is  especially  relevant  since  many  desorption  experiments  rely  on  permeation 
to  supply  hydrogen  to  the  surface.  The  research  conducted  here  has  shown  that 
the  absorption  energy  barrier  contributes  significantly  to  the  translational  energy  of 
desorbing  molecules.  Predictions  of  desorption  influenced  by  permeation  and  diffusion 
must  therefore  include  the  absorption  energy  barrier.  In  addition,  it  has  been  shown 
here  that  this  barrier  has  little  impact  on  the  final  rotational  and  vibrational  energies 
of  desorbing  molecules. 

Even  though  desorption  results  are  discussed  in  reference  to  the  H2-Cu(100)  gas- 
surface  system,  any  chemisorptive  system  with  a  Langmuir-Hinshelwood  mechanism 
could  be  applied  to  the  thermal  desorption  model.  Hence  this  work  is  universal  in 
nature  for  all  similar  systems.  The  results,  of  course,  will  depend  heavily  on  the 
accuracy  of  the  potential  energy  surface  chosen  (or  developed)  to  describe  the  gas- 
surface  interactions.  However,  with  an  accurate  potential  energy  surface,  the  physics 
of  thermal  desorption  can  be  simulated  for  any  similar  system.  Besides,  other  systems 
with  different  mechanisms  can  still  utilize  the  novel  non-dimensionalization  scheme 
presented  here,  since  the  non-dimensionalization  is  universal  for  any  classical  or  quasi- 
classical  trajectory  method. 

There  are  many  avenues  one  could  take  in  furthering  the  thermal  desorption  re¬ 
search  presented  here.  One  path  forward  would  be  to  include  adsorption  and  scat¬ 
tering  modeling  for  a  chemisorptive  system.  Since  desorption  is  the  microscopically- 
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inverse  process  of  adsorption,  it  would  not  be  difficult  to  do  so.  Initial  conditions  for 
scattering  and  adsorption  would  actually  be  simpler  to  define  than  for  desorption  be¬ 
cause  in  a  Direct  Simulation  Monte  Carlo  simulation,  one  would  know  translational 
and  internal  energies  of  impinging  molecules.  One  could  then  follow  the  molecule 
from  the  gas  phase  as  it  interacts  with  the  solid  via  the  potential  energy  surface.  If 
the  molecule  reverses  course  and  travels  back  to  the  gas  domain,  then  it  scatters. 
Otherwise,  it  adsorbs  to  the  surface.  In  this  manner,  scattering  distributions  could 
be  determined,  to  include  angular,  translational,  rotational,  and  vibrational.  In  ad¬ 
dition,  adsorption  probabilities  could  be  calculated.  Implementing  adsorption  and 
scattering  into  Direct  Simulation  Monte  Carlo  would  then  be  straightforward  and 
beneficial. 

Besides  directly  modeling  adsorption  and  scattering  for  a  chemisorptive  gas-surface 
system,  this  research  could  be  extended  by  considering  any  of  the  other  surface- 
effect  mechanisms.  These  chemical  pathways  include  Eley-Rideal,  field-induced,  ion- 
induced,  photo-induced,  and  electron-stimulated  desorption.  An  entire  suite  of  rele¬ 
vant  surface  effects  could  be  constructed  over  time  to  provide  accurate  modeling  in 
Direct  Simulation  Monte  Carlo  applications.  For  example,  in  the  field  of  laser-surface 
effects,  photo-induced  desorption  and  ablation  is  a  field  of  considerable  interest.  Per¬ 
haps  similar  methods  to  those  presented  in  this  work  could  be  developed  to  accurately 
model  dynamic  laser-surface  interactions. 

Eley-Rideal  desorption  is  another  promising  mechanism  that  could  be  investigated. 
Instead  of  modeling  two  adsorbed  atoms  on  the  surface,  one  could  model  a  single 
adatom  and  a  gas-phase  atom.  Initial  conditions  for  the  adatom  could  be  taken 
from  a  Maxwell- Boltzmann  distribution  at  the  surface  temperature.  For  the  gas- 
phase  atom,  its  initial  conditions  would  be  known  from  its  previous  interactions  and 
collisions  within  the  Direct  Simulation  Monte  Carlo  program.  One  would  need  to 


114 


determine  how  Keck’s  method  would  apply  under  these  conditions,  if  at  all.  Also,  a 
different  potential  energy  surface  would  need  to  be  created,  or  found  in  the  literature. 

Electron-stimulated  desorption  would  also  be  a  fascinating  topic  to  address.  In 
high-power  microwave  applications,  electron-stimulated  desorption  plays  a  role  in 
contaminating  the  plasma  and  inducing  pulse-shortening.  One  could  model  electron- 
stimulated  desorption  in  the  following  manner.  Take  an  impinging  electron  whose 
velocity  is  known.  The  electron  excites  a  surface  adatom  (or  admolecule)  to  some 
degree.  Then  that  adatom  travels  under  the  influence  of  the  potential  energy  surface 
to  either  desorb,  or  to  remain  on  the  surface,  albeit  at  a  new  location.  A  challenge  here 
would  be  to  incorporate  current  electron-stimulated  desorption  theory  into  a  Direct 
Simulation  Monte  Carlo  framework,  and  then  to  efficiently  model  the  desorption  of 
the  excited  adatom. 

Ablation  is  another  natural  extension  of  this  work.  Similar  to  activated  desorption, 
ablation  requires  the  overcoming  of  an  energy  barrier.  Even  though  there  would  be 
imperfections  in  a  real  surface  undergoing  ablation  (e.g.  steps)  one  could  begin  model¬ 
ing  ablation  by  assuming  that  each  ablation  event  occurs  under  the  same  conditions. 
For  example,  each  ablation  event  could  be  assumed  to  result  from  a  surface  atom 
on  a  perfectly-clean  surface  with  no  imperfections.  Ablation  would  be  induced  by 
energy  deposition,  originating  from  collisions  with  the  surface  or  conductive  heating. 
Ablation  could  therefore  be  considered,  at  least  at  first,  as  a  temperature-dependent 
process.  Initial  conditions  of  an  ablating  atom  would  then  come  from  a  Maxwell- 
Boltzmann  distribution  at  the  (local)  surface  temperature. 

In  conclusion,  this  research  both  significantly  contributes  to  the  state-of-the-art, 
and  lays  a  solid  foundation  for  future  developments  in  surface-effects  modeling  for 
Direct  Simulation  Monte  Carlo  applications.  The  effects  investigated  here  (scattering, 
adsorption,  and  Langmuir-Hinshelwood  thermal  desorption)  represent  only  a  subset 
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of  the  many  mechanisms  by  which  a  gas  interacts  with  a  solid.  Other  important 
mechanisms  that  conld  be  investigated  in  future  work  are  Eley-Rideal  desorption, 
electron-stimulated  desorption,  and  ablation.  By  modifying  and  building  upon  the 
approaches  developed  here,  the  successful  modeling  of  other  surface  effects  could  also 
be  accomplished  and  incorporated  into  Direct  Simulation  Monte  Carlo  applications. 
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Appendix  A.  List  of  Symbols 


Symbol 

a 


do 

o.\ 

0-2 

03 

04 

ai 

aQ 

A 

A 


A (2) 

Ai 

A  band 

Ac. 


Acap 


A„ 


Definition 

Average  equilibrium  length  of  the  unit  cube  [A], 
Fitting  coefficient  for  Vrep  [eV] 

Variable  in  the  generalized  quadratic  formula  [*] 
Outer  radius  [m] 

Bohr  radius  [0.529  A] 

Fitting  coefficient  for  14tt  [A  X] 

Fitting  coefficient  for  14tt  [A  2] 

Fitting  coefficient  for  14tt  [A  3] 

Fitting  coefficient  for  14tt  [A  X] 

Average  equilibrium  length  of  the  unit  cuboid  [A] 
Accommodation  coefficient  for  property  Q  [*] 
Area  of  the  surface  unit  cell  [A2] 

Unit  vector  of  a  link  [*] 

Fluxal  area  [m2] 

Average  quadrupole  alignment  parameter  [*] 

BR  potential  adjustable  parameter  [kJ  mol”1] 
Surface  area  of  a  spherical,  latitudinal  band  [m2] 
Surface  area  of  a  spherical  cap  [m2] 

BR  potential  adjustable  parameter  [kJ  mol”1] 
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Symbol  Definition 

b  Fitting  coefficient  for  Vrev  [A  4] 

Variable  in  the  generalized  quadratic  formula  [*] 
Inner  radius  [m] 

bj  aq-intercept  in  the  aio^-plane  [m] 

c  Incident  velocity  [m  s^1] 

Variable  in  the  generalized  quadratic  formula  [*] 


Co 

Fitting 

coefficient 

Cl 

Fitting 

coefficient 

Cll 

Fitting 

coefficient 

Cm 

Fitting 

coefficient 

Cnn 

Fitting 

coefficient 

Clllll 

Fitting 

coefficient 

C11112 

Fitting 

coefficient 

Cm2 

Fitting 

coefficient 

C11122 

Fitting 

coefficient 

C112 

Fitting 

coefficient 

C1122 

Fitting 

coefficient 

C11222 

Fitting 

coefficient 

C12 

Fitting 

coefficient 

C122 

Fitting 

coefficient 

C1222 

Fitting 

coefficient 

for  the  three-body  potential  [eV] 
for  the  three-body  potential  [eV  A  '] 
for  the  three-body  potential  [eV  A  "] 
for  the  three-body  potential  [eV  A  3] 
for  the  three-body  potential  [eV  A  4] 

o  _5 

for  the  three-body  potential  [eV  A 

o  _5 

for  the  three-body  potential  [eV  A 
for  the  three-body  potential  [eV  A  4] 

o  _5 

for  the  three-body  potential  [eV  A 
for  the  three-body  potential  [eV  A  3] 
for  the  three-body  potential  [eV  A  4] 

o  _5 

for  the  three-body  potential  [eV  A 
for  the  three-body  potential  [eV  A  "] 
for  the  three-body  potential  [eV  A  3] 
for  the  three-body  potential  [eV  A  4] 
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Symbol 

C12222 

c2 

C22 

C-222 

C2222 

C-22222 

C6 

C71 

d 

Ct 

4 

(4) 

C-w 

d 

do 

dA 

dE 

dE* 

du! 

dv' 


Definition 

o  _5 

Fitting  coefficient  for  the  three-body  potential  [eV  A 
Fitting  coefficient  for  the  three-body  potential  [eV  A  *] 

Fitting  coefficient  for  the  three-body  potential  [eV  A  “] 

Fitting  coefficient  for  the  three-body  potential  [eV  A  3] 

Fitting  coefficient  for  the  three-body  potential  [eV  A  4] 

o  _  5 

Fitting  coefficient  for  the  three-body  potential  [eV  A 
BR  potential  adjustable  parameter  [kJ  A6  mol-1] 

Incident  normal  velocity  [m  s-1] 

Post-collisional  normal  velocity  [m  s-1] 

Incident  tangential  velocity  [m  s-1] 

Post-collisional  tangential  velocity  [m  s-1] 

Average  post-collisional  tangential  velocity  [m  s-1] 

Wall  velocity  [m  s-1] 

Membrane  thickness  [m] 

Exponent  of  curve  fit,  cos d(6)  [*] 

Equilibrium  distance  between  surface  atoms  [A] 

Infinitesimal  area  in  the  geometric  representation  of  CLL  [*] 
Inhnitessimal  total  energy  [J] 

Inhnitessimal  (non-dimensionalized)  total  energy  [*] 

Inhnitessimal  post-collisional  tangential  velocity  component  [m  s-1] 
Inhnitessimal  post-collisional  tangential  velocity  component  [m  s-1] 
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Symbol 

dw' 

D 

D0 

D 

De 

E 

E;\(]^ 

Ed 

rp* 

■^max 

Es 

EL 


En 


E 


ELn 


ELt 


E, 

/ 


fc 


fd 

I'e 


Ie* 


Definition 

Infinitessimal  post-collisional  normal  velocity  component  [m  s-1] 
Diffusion  coefficient  [m2  s-1] 

Maximum  diffusion  coefficient  [m2  s-1] 

Step  direction  [*] 

Fitting  coefficient  for  Vatt  [eV] 

Total  energy  [J] 

Energy  of  adsorption  [kJ  mol-1] 

Activation  energy  for  diffusion  [J] 

Discrete  approximation  for  +oo  [*] 

Activation  energy  for  solubility  [J] 

Post-collisional  total  energy  [kJ  mol-1] 

Incident  translational  energy  [kJ  mol-1] 

Post-collisional  translational  energy  [kJ  mol-1] 

Post-collisional  normal  translational  energy  [kJ  mol-1] 
Post-collisional  tangential  translational  energy  [kJ  mol-1] 
Activation  energy  for  permeability  [J] 

Fugacity  [Pa] 

Switching  function  [*] 

Switching  function  [*] 

MB  energy  distribution  [J-1] 

MB  (non-dimensionalized)  energy  distribution  [*] 
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Symbol 

/max 

ft 

/ts 

ffs 

go 

9o 


9o 


9 1 
9c 

G 

G(r) 

G 

G* 

Go 

h 


It'S 

H 

Hnew 

Hold 

H* 


Definition 

Maximum  value  of  /ts  [*] 

Probability  density  of  times  between  successive  events  [*] 
Truncated  MB  (non-dimensionalizd)  energy  distribution  [*] 
/ts  normalized  by  /max  [*] 

Gradient  of  V  at  qo  [J  m"1] 

g0  at  step  i  [J  m"1] 

kth  element  of  go  [J  m-1] 

Pre-convergence  threshold  parameter  [J  m-1] 

Threshold  parameter  [J  m_1] 

Cyclic  parameter  [A  X] 

Damping  function  [*] 

Negative  gradient  of  V  at  pn  [J  m-1] 

Projected  gradient  [J  m_1] 

Component  of  g0  along  U  [J  m-1] 

Planck  constant  [6.626  x  10-34  J  s] 

Height  of  cylinder  [nr] 

BR  potential  adjustable  height  parameter  [A] 

Hessian  matrix  [J  nr"2] 

Current  value  of  H  [J  nr"2] 

Previous  value  of  H  [J  nr"2] 

Non-dimensionalized  Hessian  matrix  [*] 
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Symbol 


H' 


Ho 

Hi 

H 

H* 


[[)■ 


Hoo 

Hw 

Hu 

H BiW 


Hj,k 


TJ* 

hj,k 


H„ 


n 


i 

I 

I 

Io 

J 

J 


Definition 

Mass-weighted  Hessian  matrix  [J  kg-1  m-2] 

H  evaluated  at  qo  [J  m-2] 

H0  at  step  i  [J  m”2] 

Hamiltonian  [J] 

Non-dimensionalized  Hamiltonian  [*] 

(j,  k)  element  of  H0  [J  m-2] 

Plane-wave  diffraction  function  [A  X] 
Plane-wave  diffraction  function  [A  X] 
Plane-wave  diffraction  function  [A  X] 
Plane-wave  diffraction  function  [A  X] 

(j,  k)  element  of  H  [J  m-2] 

(j,  k)  element  of  H*  [*] 

(j,  k )  element  of  H'  [J  kg-1  m-2] 

Vibrational  term  of  the  Hamiltonian  1~L  [J] 
Reduced  Planck  constant  [1.055  x  10-34  J  s] 
Imaginary  number  a/— 1  [*] 

Reduced  moment  of  inertia  [kg  m2] 

Identity  matrix  [*] 

Modified  Bessel  Function  of  the  First  Kind  [*] 
Rotational  number  [*] 

Flux  of  molecules  [molecules  s-1  m-2] 
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Symbol 

J 

J* 

J 


Ji 


Jz 


Jz 


k 


K 

k 

h 

h 

U 

^new 

kold 

C 

m 


m  i 


m2 


mj 

M 


Definition 

Rotational  angular  momentum  [kg  m  s-1] 
Non-dimensionalized  rotational  angular  momentum  [*] 
Rotational  angular  momentum  vector  [kg  m  s^1] 
Rotational  angular  momentum  vector  of  atom  i  [kg  m  s-1] 
component  of  J  [kg  nr  s^1] 

^-component  of  J f  [kg  nr  s_1] 

Boltzmann  constant  [8.315  x  10“3  kJ  nrol-1  K_1], 
Boltzmann  constant  [1.381  x  10”23  J  K"1] 

Factor  in  the  Kisliuk  adsorption  model  [*] 

Threshold  parameter  [nr] 

Threshold  parameter  [nr] 

Threshold  parameter  [nr] 

Updated  value  of  l2  [nr] 

Previous  value  of  l2  [nr] 

Lagrangian  [J] 

Molecular  mass  [Mg  mol-1], 

Mass  of  an  atom  [kg] 

Mass  of  atom  1  [kg] 

Mass  of  atom  2  [kg] 

Slope  of  a  line  [*] 

Total  mass  [kg] 
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Symbol 

M1 

M2 

W-m 

nv 

N 


N colls 

p* 


p'*. 

■L  C.l 


Ph 

Ph-i 

Ph+i 

Ph new 

Phm 

Pi 

Pi 

Pi 


Pr 

Pr 

Pr 

Pr 


Definition 

Starting  configuration 

Ending  configuration 

Exponent  of  the  modified  energy  [*] 

Vibrational  number  [*] 

Number  of  atoms  in  the  molecule  [atoms] 

Number  of  adatoms  [adatoms] 

Number  of  collisions  [*] 

New  point  for  replacing  p h 
Non-dimensionalized  MWC  momenta  [*] 

Approximate  location  of  TS 
Path  point  adjacent  to  pn 
Path  point  adjacent  to  pn 
Current  value  of  pn 
Previous  value  of  pn 
Generalized  momentum  [kg  m  s-1] 

Generalized  force  [kg  m  s-2] 

Momentum  vector  of  atom  i  [kg  m  s'1] 

Momentum  in  the  r-direction  [kg  m  s_1] 
Non-dimensionalized  momentum  in  the  r-direction  [*] 
Force  in  the  r-direction  [kg  m  s-2] 
Non-dimensionalized  force  in  the  r-direction  [*] 
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Symbol 


P 


U,l 


Px 

P*x 


Px  1 


V 


* 

XI 


p 


x\ 


p 


/  * 
Xl 


Px  2 


P  X2 
P'x2 

p'* 

rx2 

Px 

P*x 


p. 


X,l 


Py 


P 


* 

y 


Py  i 
P*y  i 


p 


yi 


P 


y  i 


Py2 


Definition 

Non-dimensionalized  momentum  in  the  MWC  Hessian  eigenspace  [*] 
Momentum  in  the  x-direction  [kg  m  s-1] 

Non-dimensionalized  momentum  in  the  x- direction  [*] 

Momentum  in  the  Xi -direction  [kg  m  s”1] 

Non-dimensionalized  momentum  in  the  xi-direction  [*] 
Mass-weighted  momentum  in  the  x  i  -direction  [kg1/2  m  s^1] 
Non-dimensionalized  value  of  p1  [*] 

Momentum  in  the  ^-direction  [kg  m  s”1] 

Non-dimensionalized  momentum  in  the  ^-direction  [*] 
Mass-weighted  momentum  in  the  x'2-direction  [kg1/2  m  s^1] 
Non-dimensionalized  value  of  p'X2  [*] 

Force  in  the  x- direction  [kg  m  s-2] 

Non-dimensionalized  force  in  the  x- direction  [*] 
x-component  of  Pi  [kg  m  s~x] 

Momentum  in  the  ^/-direction  [kg  m  s_1] 

Non-dimensionalized  momentum  in  the  y-direction  [*] 

Momentum  in  the  yi-direction  [kg  m  s^1] 

Non-dimensionalized  momentum  in  the  yi-direction  [*] 
Mass-weighted  momentum  in  the  y!-direction  [kg1/2  m  s”1] 
Non-dimensionalized  value  of  p’  [*] 

Momentum  in  the  y2-direction  [kg  m  s~x] 
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Symbol 


Definition 


P 


V2 


P 


V2 


P 


V2 


Py 


P 


* 

y 


Py,i 


Pz 


P 


* 

2: 


Pzi 


P*Z1 
P'z ! 


V 


/  * 
Z1 


Pz2 

P*z2 


p 


Z  2 


p 


z  2 


Pz 


Pe 

Pi 

Pe 


Non-dinrensionalized  momentum  in  the  ^-direction  [*] 
Mass-weighted  momentum  in  the  ^-direction  [kg1/2  m  s”1] 
Non-dimensionalized  value  of  p'y2  [*] 

Force  in  the  ^-direction  [kg  m  s-2] 

Non-dimensionalized  force  in  the  y-direction  [*] 
^-component  of  Pi  [kg  nr  s^1] 

Momentum  in  the  ^-direction  [kg  m  s_1] 
Non-dimensionalized  momentum  in  the  ^-direction  [*] 
Momentum  in  the  Z\ -direction  [kg  nr  s^1] 
Non-dimensionalized  momentum  in  the  zi -direction  [*] 
Mass-weighted  momentum  in  the  Zi -direction  [kg1/2  nr  s^1] 
Non-dimensionalized  value  of  p'zi  [*] 

Momentum  in  the  ^-direction  [kg  nr  s^1] 
Non-dimensionalized  momentum  in  the  ^-direction  [*] 
Mass-weighted  momentum  in  the  ^-direction  [kg1/2  nr  s~x] 
Non-dimensionalized  value  of  p'  [*] 

Force  in  the  ^-direction  [kg  nr  s-2] 

Non-dimensionalized  force  in  the  ^-direction  [*] 

Momentum  in  the  ^-direction  [kg  nr  s_1] 
Non-dimensionalized  momentum  in  the  0-direction  [*] 
Force  in  the  0-direction  [kg  nr  s-2] 
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Symbol 

Definition 

P*e 

Non-dimensionalized  force  in  the  ^-direction  [*] 

P<p 

Momentum  in  the  ^-direction  [kg  m  s-1] 

Pi 

Non-dimensionalized  momentum  in  the  ^-direction  [*] 

t>4> 

Force  in  the  0-direction  [kg  m  s-2] 

Pi 

Non-dimensionalized  force  in  the  0-direction  [*] 

P 

Gas  partial  pressure  [Pa] 

P(si,s2) 

Least-squares  fit  for  the  three-body  potential  [eV] 

P(u  — >  u') 

Scattering  kernel  [*] 

P(v  — >  v') 

Scattering  kernel  [*] 

P(w  — y  w') 

Scattering  kernel  [*] 

Qo 

Initial  configuration  vector  in  Cartesian  coordinates  [m 

<7o 

qo  at  step  i  [m] 

Qt 

Generalized  coordinate  [m] 

Qi 

Generalized  velocity  [m  s"1] 

Qi 

Generalized  acceleration  [m  s~2] 

Qm 

Factor  in  the  MK  adsorption  model  [*] 

Q 

Molecular  property  in  question 

Total  permeation  rate  [molecules  s-1] 

Q 

Component  of  the  step  A q  along  U  [m] 
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Symbol 

Definition 

r 

Intermolecular  distance  [A]  [m] 

Desorption  rate  [molecules  s_1] 

Detector  radius  [m] 

r 

Non-dimensionalized  r-coordinate  [*] 

r 

Unit  vector  in  the  r-direction  [*] 

r 

Time-rate-of-change  of  the  r-coordinate  [m  s-1] 

y* 

Non-dimensionalized  time-rate-of-change  of  the  r-coordinate  [* 

r0 

Equilibrium  intermolecular  distance  [A], 

Fitting  coefficient  for  the  three-body  potential  [A] 

T\ 

BR  potential  adjustable  parameter  [A] 

r3  b 

Fitting  coefficient  for  the  three-body  potential  [A] 

rc 

Distance  in  geometric  representation  of  CLL  [*] 

re 

Fitting  parameter  for  Vatt  [A] 

rn  s 

Radius  of  the  hypersphere  [m] 

Vi 

Distance  between  gas  molecule  and  ith  bulk  molecule  [A] 

K 

Adjusted  value  of  rt  [A] 

Vi 

Radius  vector  of  atom  i  [m] 

TIef 

Reaction  zone  parameter  [A] 

rref,2 

Reaction  zone  parameter  [A] 

T  x,i 

^-component  of  ry  [m] 

rv,i 

^/-component  of  ry  [m] 

128 


Symbol 


Definition 


R 


n 

Ri 

R'2 

r3 

R4 

R-c 

RcS 

Ri 

Rsp 

S 


Sl 


s  2 


So 


S(0) 

t 

(t) 

t * 


Uniform  random  number  between  0  and  1  [*] 
Distance  from  the  detector  to  the  surface  [m] 
Radius  parameter  [A] 

Uniform  random  number  between  0  and  1  [*] 
Uniform  random  number  between  0  and  1  [*] 
Uniform  random  number  between  0  and  1  [*] 
Uniform  random  number  between  0  and  1  [*] 
Critical  effective  interaction  radius  [A] 
Effective  interaction  radius  [A] 

Uniform  random  number  between  0  and  1  [*] 
Specific  gas  constant  [J  kg"1  K”1] 

Solubility  [molecules  m-3  Pa -1/2] 

Variable  for  the  three-body  potential  [A] 
Variable  for  the  three-body  potential  [A] 
Initial  adsorption  probability  [*] 

Maximum  solubility  [molecules  m-3  Pa-1/2] 
Adsorption  probability  onto  the  adlayer  [*] 
Adsorption  probability  [*] 

Time  [s] 

Average  time  between  events  [s] 
Non-dimensionalized  time  [*] 


129 


Symbol 

tmp 

T 


T 

* j — * 

T 

T 


U 

u' 


Ump 

U 


U 


v 

V 

V ' 

V 
V* 
V0 

Voooo 

^0010 


Definition 

Most-probable  travel  time  [s] 

Pseudo-tangent  to  the  path 
Kinetic  energy  [J] 

Non-dimensionalized  kinetic  energy  [*] 

Critical  temperature  [K] 

Surface  temperature  [K] 

Triple-point  temperature  [K] 

A  tangential  velocity  component  [m  s-1] 

A  post-collisional  tangential  velocity  component  [m  s"1] 
jth  Hessian  eigenvector  [*] 

Most-probable  speed  [m  s-1] 

Vector  in  calculating  the  BFGS  algorithm  [m  J^1] 
Unitary  matrix  that  diagonalizes  H  [*] 

A  tangential  velocity  component  [m  s-1] 

Velocity  [m  s_1] 

A  post-collisional  tangential  velocity  component  [m  s^1] 
Potential  energy  [J] 

Non-dimensionalized  potential  energy  [*] 

Potential  energy  at  some  initial  configuration  [J] 

6D  expansion  coefficient  [eV  A] 

6D  expansion  coefficient  [eV  A] 
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Symbol 

boon 

Voob 

VoOh 

y^t 


V2 


VA 

v2 


V* 


V2000 

V2010 

V2011 

^2el0 

V20  b 
V20  h 
ViOt 

^2e6 


^3 

tM 

K3 

T/S 

K3 

b^D 

baa 


Definition 

6D  expansion  coefficient  [eV  A] 

2D  expansion  coefficient  [eV] 

2D  expansion  coefficient  [eV] 

2D  expansion  coefficient  [eV] 

Two-body  potential  [eV] 

Two-body  potential  for  H2  as  it  desorbs  [eV] 
Two-body  potential  for  two  H  atoms  [eV] 

6D  expansion  coefficient  [eV  A] 

6D  expansion  coefficient  [eV  A] 

6D  expansion  coefficient  [eV  A] 

6D  expansion  coefficient  [eV  A] 

2D  expansion  coefficient  [eV] 

2D  expansion  coefficient  [eV] 

2D  expansion  coefficient  [eV] 

2D  expansion  coefficient  [eV] 

Three-body  potential  [eV] 

Three-body  potential  for  H2  as  it  desorbs  [eV] 
Three-body  potential  for  two  H  atoms  [eV] 
Wiesenekker  6D  PES  [eV] 

Two-body  attractive  potential  [eV] 

2D  PES  [eV] 
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Symbol 
Vbh90 
Vbt9  0 

Hil40 

Vrr 

Vhb90 

VhbUO 

Vb 

Rep 

Vtb90 

Rbl40 

V(z: ave) 

w 

w' 

W 

X 

X* 

X 

X* 

X\ 


Definition 
2D  PES  [eV] 

2D  PES  [eV] 

2D  PES  [eV] 

BR  potential  [kJ  mol-1] 

2D  PES  [eV] 

2D  PES  [eV] 

LJ  potential  [kJ  mol-1] 

Morse  potential  [kJ  mol-1] 

Two-body  repulsive  potential  [eV] 

2D  PES  [eV] 

2D  PES  [eV] 

BR  potential  term  [kJ  mol-1] 

Normal  velocity  component  [m  s-1] 

Post-collisional  normal  velocity  component  [m  s-1] 

BR  potential  adjustable  parameter  [kJ  mol-1] 

^-coordinate  (parallel  to  the  surface)  [m] 

Reaction  order  for  desorption  [*] 

Non-dimensionalized  re-coordinate  [*] 

Time-rate-of- change  of  the  x-coordinate  [m  s-1] 
Non-dimensionalized  time-rate-of-change  of  the  x-coordinate  [*] 
x-coordinate  of  atom  1  [m] 
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Symbol 

X\ 

x2 

X-2 

Xg 

Xi 


y 

y* 

y 

y* 

yi 

ij\ 

y2 

h 

y3 


yi 


*?(M) 

Yim{OA) 


z 


* 


Definition 

Time-rate-of- change  of  the  x-coordinate  of  atom  1  [m  s_1] 
x-coordinate  of  atom  2  [m] 

Time-rate-of- change  of  the  x-coordinate  of  atom  2  [m  s”1] 

A  coordinate  of  the  gas  molecule  [A] 

A  coordinate  of  the  surface  molecule  [A] 
x-coordinate  of  atom  i  [m] 

//-coordinate  (parallel  to  the  surface)  [m] 

Non-dimensionalized  ^/-coordinate  [*] 

Time-rate-of- change  of  the  //-coordinate  [m  s-1] 
Non-dimensionalized  time-rate-of-change  of  the  //-coordinate  [*] 
//-coordinate  of  atom  1  [m] 

Time-rate-of-change  of  the  //-coordinate  of  atom  1  [m  s”1] 
//-coordinate  of  atom  2  [m] 

Time-rate-of-change  of  the  //-coordinate  of  atom  2  [m  s”1] 

A  coordinate  of  the  gas  molecule  [A] 

A  coordinate  of  the  surface  molecule  [A] 

//-coordinate  of  atom  i  [m] 

Modified  spherical  harmonic  [*] 

Spherical  harmonics  of  degree  l  and  order  m  [*] 

Height  above  the  surface  [A]  [m] 

Non-dimensionalized  ^-coordinate  [*] 
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Symbol 

z 

z * 

z0 

Zi 

z-1 

Z2 

Z2 

Z3b 

Ze 

Z9 

„ave 

"g 

~ave 

Zgs 

Zi 

Aef 

Zref,2 


Definition 

Time-rate-of- change  of  the  ^-coordinate  [m  s-1] 
Non-dimensionalized  time-rate-of-change  of  the  ^-coordinate  [*] 
Fitting  coefficient  for  the  three-body  potential  [A] 

Height  of  atom  1  relative  to  the  surface  [A]  [m] 
Time-rate-of-change  of  the  ^-coordinate  of  atom  1  [m  s-1] 
Height  of  atom  2  relative  to  the  surface  [A]  [m] 
Time-rate-of-change  of  the  ^-coordinate  of  atom  2  [m  s”1] 
Fitting  coefficient  for  the  three-body  potential  [A] 

Fitting  parameter  for  Vatt  [A] 

A  coordinate  of  the  gas  molecule  [A] 

Height  above  the  “local-average”  surface  [A] 

Height  of  the  “local-average”  surface  [A] 

A  coordinate  of  the  surface  molecule  [A] 

Reaction  zone  parameter  [A] 

Reaction  zone  parameter  [A] 


Greek  Symbol  Definition 

a  BR  potential  adjustable  parameter  [A  '] 

BR  potential  adjustable  parameter  [A  '] 

A  Cartesian  coordinate  [m] 

alc  ai-coordinate  of  the  center  of  the  hypersphere  [m] 
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Greek  Symbol  Definition 


CXjc 

aj,i 

1 

^n,  ad 
OLt 

^t,ad 

A 

A 

71 

72 
5 

Ag 

A9j 

A&- 

At 

Ay 

Ay* 

Ax 


A  Cartesian  coordinate  [m] 

a, -coordinate  of  the  center  of  the  hypersphere  [m] 
ay  taken  from  [m] 
ay  taken  from  Pi+i  [m] 

Accommodation  coefficient,  normal  kinetic  energy  [*] 

Adlayer  accommodation  coefficient,  normal  kinetic  energy  [*] 
Accommodation  coefficient,  tangential  kinetic  energy  [*] 

Adlayer  accommodation  coefficient,  tangential  kinetic  energy  [*] 
Generalized  mass-weighted  coordinate  [m  kg1/2] 

Generalized  mass-weighted  acceleration  [m  s-2  kg1/2] 

Fitting  coefficient  for  the  three-body  potential  [A  X] 

Fitting  coefficient  for  the  three-body  potential  [A  *] 

BR  potential  adjustable  parameter  [A  "] 

Configuration  perturbation  vector  in  Cartesian  coordinates  [m] 
A  q  with  all  elements  zero  except  for  the  jth  element  [m] 
jth  element  of  A  q[  m] 

Desorption  event  time  step  [s] 

Change  in  V  between  adsorbed  and  TS  configurations  [J] 
Non-dimensionalized  value  of  Ay  [*] 

Difference  between  x\  (£2)  and  x  [m] 

Step  taken  in  evaluating  a  central  difference  derivative  [m] 
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Greek  Symbol  Definition 


Ay 

Az 

AC 
a  e 
AC 

^LJ 

£1 

e2 

c 

Co 

v 

9 


e 

9 

9* 

0C 


Difference  between  y\  (y2)  and  y  [m] 

Difference  between  z\  (z2)  and  z  [m] 

Reaction  zone  parameter  [radians] 

Angular  sweep  of  Acap  [radians] 

Reaction  zone  parameter  [radians] 

LJ  well-depth  parameter  [kJ  mol-1] 

Morse  well-depth  parameter  [kJ  mol-1] 

Threshold  parameter  [m] 

Threshold  parameter  [*] 

Reaction  zone  parameter  [radians] 

Reaction  zone  parameter  [radians] 

Threshold  parameter  [radians] 

Adlayer  coverage  [ML], 

Elevation  angle  [radians] 

Angle  between  G  and  T  [radians] 

Desorption  angle,  measured  from  the  surface  normal  [radians] 
Non-dimensionalized  ^-coordinate  [radians] 

Unit  vector  in  the  ^-direction  [*] 

Time-rate-of- change  of  the  0-coordinate  [radians  s-1] 
Non-dimensionalized  value  of  9  [radians] 

Angle  in  geometric  representation  of  CLL  [*] 
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Greek  Symbol  Definition 


Of 

Oi 

A 

Ao 

I1 


e 


p&tt 


Prep 


CTlj 

&M 

® n 


<?t 


4> 

4> 


In-plane  scattering  angle  with  respect  to  surface  normal  [°] 
Angle  of  incidence  with  respect  to  the  surface  normal  [°] 
Optimal  step  length  [m] 

Diagonal  matrix  of  the  eigenvalues  of  Ho  [*] 

Reduced  mass  [kg] 

Molecular  vibrational  frequency  [cycles  s_1] 

Frequency  factor  [molecules-^-1'  g-1] 

Reaction  zone  parameter  [radians] 

Angle  between  the  z-axis  and  J  [radians] 

Reaction  zone  parameter  [radians] 

Variable  for  14tt  [A] 

Variable  for  Vrep  [A] 

LJ  zero-potential  intermolecular  distance  [A] 

Morse  potential  parameter  [A  *] 

Accommodation  coefficient,  normal  momentum  [*] 
Accommodation  coefficient,  tangential  momentum  [*] 
Azimuthal  angle  [radians] 

Concentration  of  molecules  [molecules  m-3] 
Non-dimensionalized  0-coordinate  [radians] 

Unit  vector  in  the  ^-direction  [*] 

Time-rate-of- change  of  the  0-coordinate  [radians  s^1] 
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Greek  Symbol 

</>* 

<f> 

$0 

<!>(?) 

l 

(bQ 

r 

cbQ 

W 

X 

X2 

o 


Definition 

Non-dimensionalized  value  of  [radians] 
Permeability  [molecules  m-1  s_1  Pa _1/2] 

Maximum  permeability  [molecules  m"1  s_1  Pa _1/2] 
BR  potential  weighting  function  [*] 

Incident  flux  of  Q  [cm-2  s"1] 

Reflected  flux  of  Q  [cm-2  s-1] 

Fully-accommodated  reflected  flux  of  Q  [cm"2  s_1] 
Reaction  zone  parameter  [radians] 

Reaction  zone  parameter  [radians] 

Solid  angle  [steradians] 
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Appendix  B.  List  of  Acronyms 


Acronym 

Definition 

ID 

One-Dimensional 

2D 

Two-Dimensional 

3D 

Three-Dimensional 

6D 

Six-Dimensional 

AFOSR 

Air  Force  Office  of  Scientific  Research 

AFRL 

Air  Force  Research  Laboratory 

BFGS 

Br  oy  den-  Flet  cher-  Goldfarb-  S  hanno 

BR 

Barker-Rettner 

CDF 

Cumulative  Distribution  Function 

CFD 

Computational  Fluid  Dynamics 

CLL 

Cercignani-Lampis-Lord 

CT 

Classical  Trajectory 

DFT 

Density  Functional  Theory 

DIET 

Desorption  Induced  by  Electronic  Transitions 

DSMC 

Direct  Simulation  Monte  Carlo 

EM 

Electromagnetic 

ER 

Eley-R  ideal 

ESD 

Electron-Stimulated  Desorption 

FCC 

Face- Centered  Cubic 

FID 

Field-Induced  Desorption 
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Acronym 

Definition 

GGA 

Generalized  Gradient  Approximation 

GW 

Gaussian  Weighting 

HPM 

High-Power  Microwave 

ICEPIC 

Improved  Concurrent  Electromagnetic  Particle-In-Cell 

IID 

Ion-Induced  Desorption 

LEPS 

London-Eyring-Polyani-Sato 

LH 

Langmuir- Hinshelwood 

LJ 

Lennard-Jones 

MB 

Maxwell- Boltzmann 

MC 

Monte  Carlo 

MD 

Molecular  Dynamics 

MEMS 

Micro- Electro-Mechanical  Systems 

MK 

Modified  Kisliuk 

MKS 

Modified  Kisliuk  with  Scattering 

ML 

Monolayer 

MS 

Multi-Stage 

MWC 

Mass- Weighted  Cartesian 

PES 

Potential  Energy  Surface 

PID 

Photo- Induced  Desorption 

QCT 

Quasi-Classical  Trajectory 

RS 

Rigid  Surface 
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Acronym 

Definition 

SEE 

Secondary  Electron  Emission 

sws 

Slow  Wave  Structure 

TPD 

Temperature-Programmed  Desorption 

TS 

Transition  State 

TST 

Transition  State  Theory 

TTPD 

Threshold  Temperature-Programmed  Desorption 

UHV 

Ultra-High  Vacuum 

XHV 

Extreme-High  Vacuum 

ZPE 

Zero-Point  Energy 
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